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ABSTRACT:  This  report  documents  the  development  of  a  verified  2-D  numerical  model  of  the  hydro¬ 
dynamics  and  salinity  of  Biscayne  Bay,  FL.  The  computer  code  employed  for  this  study  was  TABS-MDS 
(multidimensional  sediment)  (formerly  called  RMA10-WES).  The  model  was  calibrated  and  verified 
against  an  extensive  set  of  hydrodynamic  and  salinity  field  data  collected  in  a  previous  effort  performed 
by  the  U.S.  Army  Engineer  Research  and  Development  Center,  Coastal  and  Hydraulics  Laboratory,  in 
conjunction  with  the  Biscayne  National  Park.  This  report  discusses  general  characteristics  of  the 
Biscayne  Bay  system,  model  development,  model  boundary  condition,  and  calibration  and  verification 
results. 
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The  model  development  investigation  reported  herein  was  conducted  by  the 
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Director,  Information  Technology  Laboratory  and  Dr.  William  H.  McAnally, 
retired,  CHL,  were  also  very  valuable. 

COL  James  R.  Rowan,  EN,  was  Commander  and  Executive  Director  of 
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Introduction 


Background 

Biscayne  Bay,  near  Miami,  FL,  is  a  popular  tourist  destination,  an  important 
commercial  fishing  resource,  and  a  valuable  habitat  for  several  endangered 
species  of  flora  and  fauna.  Therefore,  the  U.  S.  Army  Engineer  District, 
Jacksonville,  as  well  as  local  sponsors  and  partners  such  as  the  Miami-Dade 
Department  of  Environmental  Resource  Management  (DERM),  Biscayne 
National  Park  (BNP),  and  the  South  Florida  Water  Management  District 
(SFWMD)  are  interested  in  developing  data  sets  and  numerical  models  that  can 
aid  in  the  study  and  management  of  Biscayne  Bay  circulation,  salinity,  and  water 
quality.  For  this  purpose,  the  U.  S.  Army  Engineer  Research  and  Development 
Center  (ERDC),  Coastal  and  Hydraulics  Laboratory  (CHL),  was  tasked  with 
developing  a  physics  based  numerical  model  of  Biscayne  Bay. 

The  study  of  circulation,  salinity,  and  water  quality  patterns  in  a  system  such 
as  Biscayne  Bay  is  a  complex  issue.  Physical  processes  that  impact  the  water 
quality  within  the  system  vary  both  spatially  and  temporally .  F actors  that  deter¬ 
mine  circulation  and  salinity  patterns  in  the  system'include:  the  bathymetry  and 
geometry  of  the  navigation  channels,  interconnecting  canals  and  inlets;  astro¬ 
nomical  tide-induced  currents;  wind-induced  currents;  freshwater  inflow;  rainfall 
and  evaporation;  and  vertical  density  gradients  resulting  from  salinity  stratifica¬ 
tion.  This  numerical  model  study  was  proposed  and  undertaken  primarily  to 
develop  numerical  modeling  tools  to  aid  in  the  further  assessment  of  the  impact 
of  these  and  other  features  on  the  system.  Particular  emphasis  was  placed  on  the 
use  of  the  tools  to  assess  the  impact  of  changing  freshwater  inflows  on  the  hydro¬ 
dynamics  and  salinity  of  the  bay.  The  companion  field  data  collection  effort  is 
discussed  by  McAdory  et  al.  (2002). 


Location  and  description  of  system 

Biscayne  Bay  is  a  shallow,  subtropical  marine  lagoon,  with  substantial 
estuarine  characteristics  along  the  western  shoreline.  It  is  located  along  the 
southeast  coast  of  Florida.  Biscayne  Bay  is  about  60  miles1  long  and  a  maximum 
of  about  8  miles  wide.  The  bay  extends  from  Jewfish  Creek  in  Barnes  Sound  in 
the  south  to  Dumfoundling  Bay  in  the  north.  The  Florida  mainland  bounds  the 


1  Units  of  measurement  cited  in  this  report  are  in  non- SI  units.  A  table  of  factors  for  converting 
non-SI  to  SI  units  of  measurement  is  presented  on  p.  vi. 
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bay  on  the  west,  and  the  eastern  boundaiy  is  formed  by  a  string  of,  primarily, 
coral  islands  and  shallow,  vegetated  mud  banks  from  Key  Largo  in  the  south  to 
Miami  Beach  in  the  north.  A  location  map  of  Biscayne  Bay  is  shown  in  Figure  1. 
Figure  2  gives  the  locations  of  various  significant  land  and  water  features  within 
the  bay.  Fresh  water  enters  the  bay  primarily  through  16  man-made  coastal 
salinity  control  structures,  which  allow  a  controlled  release  of  surface-water  run¬ 
off  via  a  system  of  drainage  canals  and  levees.  Additional  fresh  water  is  con¬ 
tributed  directly  through  rainfall  and  groundwater  flows. 

In  terms  of  hydrodynamics,  the  bay  can  be  divided  into  three  distinct  regions: 
the  northern,  middle,  and  southern  reaches.  The  northern  reach  (North  Bay) 
extends  northward  from  Rickenbacker  Causeway,  a  bridge  linking  Key  Biscayne 
with  the  city  of  Miami,  to  Dumfoundling  Bay.  Man-made  channels,  dredged 
areas,  and  man-made  islands  characterize  this  northern  portion  of  the  bay.  The 
circulation  in  this  area  is  not  affected  significantly  by  the  wind.  Instead,  the 
various  channels  and  canals  guide  the  tidal  exchange  with  the  Florida  Straits 
through  several  inlets.  Government  Cut  on  the  south  and  Bakers  Haulover  Inlet 
on  the  north  are  two  of  the  main  direct  connections  of  North  Bay  with  the  off¬ 
shore  area.  The  area  in  North  Bay  between  Government  Cut  and  Bakers 
Haulover  Inlet  is  characterized  by  the  existence  of  a  tidal  node  (or  standing 
wave)  where  little  horizontal  excursion  due  to  tidal  forcing  is  observed. 

The  central  reach  of  the  bay  (Central  Bay)  extends  from  Card  Sound  north¬ 
ward  to  Rickenbacker  Causeway.  The  central  portion  of  Biscayne  Bay  is 
exposed  to  direct  tidal  exchange  with  the  Florida  Straits  through  the  Safety 
Valve,  a  natural,  shallow  inlet  approximately  8  miles  wide  and  bounded  by  Key 
Biscayne  on  the  north  and  the  Ragged  Keys  and  Elliott  Key  on  the  south.  Central 
Bay  is  characterized  by  large  expanses  of  shallow,  open  water.  Consequently, 
the  wind  speed  and  direction,  along  with  the  tide,  have  a  strong  effect  on  the 
currents  in  this  region. 

The  southern  reach  of  Biscayne  Bay  (South  Bay)  is  formed  by  Card  Sound, 
Barnes  Sound,  and  adjacent  waters,  which  experience  direct  tidal  exchange  with 
the  Florida  Straits  through  a  series  of  natural  inlets  into  Card  Sound.  These  inlets 
include  Angelfish  Creek,  Broad  Creek,  and  Caesar  Creek  (known  as  the  ABC’s). 
Barnes  Sound  is  isolated  from  direct  ocean  tidal  influence.  This  southern  portion 
of  the  bay  exhibits  a  complicated  pattern  of  tidal  exchange  and  wind-driven  cur¬ 
rents.  Water  can  enter  and  leave  this  area  through  any  of  the  several  mentioned 
natural  inlets  exchanging  with  the  Florida  Straits  to  the  east,  to  the  Central  Bay  to 
the  north,  or  through  Jewfish  Creek  to  Florida  Bay  to  the  south.  Several  small 
connections  from  Manatee  Bay  to  Long  Sound  in  Florida  Bay  also  exist  through 
culverts  under  U.S.  Highway  1. 

The  present  day  salinity  of  Biscayne  Bay  is  governed  by  the  addition  of  small 
volumes,  compared  to  the  tidal  prism,  of  fresh  water  from  various  sources  (pri¬ 
marily  canals,  rainfall,  and  groundwater)  combined  with  the  influence  of  tidal 
exchange  with  the  Florida  Straits  through  various  inlets  and  channels  and  wind- 
driven  circulation.  In  general,  the  isohalines  tend  to  align  themselves  in  a  north- 
south  direction,  with  the  fresher  water  remaining  closer  to  the  western  shoreline. 
The  north-south  alignment  changes  with  inflow  patterns,  usually  seasonally,  with 
the  isohalines  shifting  seaward  or  landward,  depending  on  the  amount  of  fresh 
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water  present.  Local  conditions  of  rainfall  and  groundwater  additions,  as  well  as 
the  occurrence  of  extreme  meteorological  events  (such  as  the  landfall  of  a  major 
hurricane)  also  affect  the  distribution  of  salinity  in  the  bay. 

In  the  past,  the  bay  proper  was  fresher  than  it  is  today,  with  the  higher  iso¬ 
halines  farther  from  the  western  boundary  than  is  now  the  case.  However, 
channelization  and  impoundment  of  overland  freshwater  flows,  together  with  the 
creation  of  Bakers  Haulover  Inlet  and  the  expansion  of  Government  Cut,  have 
resulted  in  a  significant  increase  in  the  average  salinity  of  the  system.  Bellmund 
et  al.  (1999)  and  Cantillo  et  al.  (2000)  provide  good  discussions  of  the  bay. 

Although  there  is  evidence  that  Biscayne  Bay  experiences  local  intermittent 
stratification,  particularly  in  the  vicinity  of  the  canal  mouths,  the  system  is 
believed  to  be  largely  well  mixed,  primarily  the  result  of  the  influence  of  wind 
mixing  (Alleman  et  al.  1 995). 


Study  objectives 

The  goals  of  this  effort  were  to  develop  a  two-dimensional  (2-D)  finite 
element  (FE)  hydrodynamic  and  salinity  transport  model  and  to  use  the  model  to 
characterize  the  response  of  the  bay  system  salinity  and  circulation  to  various 
freshwater  discharge  scenarios.  This  report  documents  efforts  toward  the 
development  of  this  verified  2-D  numerical  model.  Results  of  the  scenario 
experiments  using  the  model  will  be  reported  in  a  later  publication.  The  hydro- 
dynamic  and  salinity  results  derived  from  the  use  of  this  verified  numerical 
model  can  be  used,  together  with  additional  water  quality  models,  to  simulate  the 
environmental  impacts  of  various  proposed  changes  to  the  regulation  of  fresh¬ 
water  discharge  to  the  system. 


Previous  work 

This  effort  is  not  the  first  attempt  at  developing  a  2-D  computational  hydro- 
dynamic  and  salinity  transport  model  of  Biscayne  Bay.  Although  several  earlier 
models  were  developed,  the  report  by  Fatt  and  Wang  (1987)  stands  as  a  major 
source  of  information  about  other  previous  models,  as  well  as  their  own. 

Fatt  and  Wang  (1987)  describe  the  results  of  using  an  FE  model  to  simulate 
salinity  transport  in  Biscayne  Bay.  Results  from  an  earlier  FE  hydrodynamic 
model,  Wang  and  Connor  (1975),  were  used  to  drive  this  salinity  transport 
model.  The  Wang  and  Connor  hydrodynamic  model  was  verified  and  docu¬ 
mented  by  Wang  (1978).  Fatt  and  Wang’s  salinity  transport  results  were  then 
compared  against  synoptic  field  measurements  collected  by  the  National  Park 
Service  (Fatt  and  Wang  1987).  These  data  were  collected  from  August  through 
October  1985.  Of  interest  are  not  only  the  results  of  the  salinity  transport 
modeling,  but  the  analysis  of  the  hydrodynamic  and  constituent  transport 
characteristics  of  the  bay  in  general. 
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The  FE  mesh  used  by  Fatt  and  Wang  for  constituent  transport  was  identical 
to  that  on  which  the  hydrodynamic  results  were  computed.  It  has  quadratic, 
triangular  elements  and  covers  the  area  from  Rickenbacker  Causeway  south  to 
the  Card  Bank.  The  western  boundary  is  land  while  the  eastern  boundary  is 
along  the  various  inlets,  including  Safety  Valve,  Bear  Cut,  Caesar’s  Creek,  etc. 

A  total  of  287  elements  were  used.  The  algorithm  by  which  the  salinity  values 
were  computed  involved  a  hybrid  Eulerian-Lagrangian  scheme  whereby  the 
transport  equation  was  decomposed  into  convective  and  diffusive  equations, 
solved  separately,  and  the  results  summed. 

The  boundary  conditions  prescribed  to  the  Fatt  and  Wang  salinity  transport 
model  were  those  at  the  canal  mouths  and  along  the  inlets  exchanging  with  the 
Florida  Straits.  The  salinity  value  assigned  to  nodes  along  ocean  inlet  boundary 
locations  was  chosen  based  on  salinity  data  gathered  at  Adams  Key  and 
Markers  1/1A.  Salinity  values  at  the  canal  mouth  locations  were  assigned  using  a 
reduction  factor,  which  compensates  for  the  mixing  of  fresh  and  salt  water  that 
occurs  in  the  canals  upstream  of  the  bay.  A  separate  reduction  factor  was  esti¬ 
mated  for  each  canal  mouth,  using  a  trial-and-error  calibration  procedure.  This 
procedure  has  the  advantage  of  compensating  for  any  vertical  mixing  which  may 
occur  in  the  canals,  thereby  mitigating  some  of  the  impact  of  using  a  vertically 
averaged  model  to  simulate  locally  stratified  flow.  However,  the  procedure 
suffers  from  being  calibrated  to  a  narrow  range  of  discharge  rates  and  durations, 
and  therefore  each  reduction  factor  is  only  applicable  when  the  canals  are  flow¬ 
ing  at  similar  discharges  and  under  similar  circumstances  to  the  ones  for  which 
the  reduction  factors  were  calculated  (Fatt  and  Wang  1987).  Rainfall  data  were 
used  as  a  model  input,  but  evaporation  and  groundwater  flows  were  considered 
of  little  consequence  and  therefore  not  included  in  the  Fatt  and  Wang  model. 

The  salinity  transport  model  was  run  for  the  period  October  20  through  25, 
1985,  using  selected  hydrodynamic  solution  fields  to  drive  it.  Comparisons  were 
made  between  model  simulations  and  field  measurements  with  contour  plots  of 
salinity.  Overall,  the  model  was  able  to  predict  the  same  general  pattern  of 
salinity  distribution  as  shown  by  the  field  measurements.  The  salinity  contours 
are  aligned  primarily  in  a  north-south  direction,  with  some  local  closure  around 
the  canal  mouths  along  the  western  shoreline. 

Also  of  interest  in  the  Fatt  and  Wang  (1987)  study  are  the  field  data.  Their 
analysis  of  the  synoptic  salinity  measurements,  at  some  48  stations  throughout 
the  central  bay,  reveals  some  interesting  facts.  Most  readings  were  taken  at  low 
or  high  tide,  at  three  locations  in  the  water  column  (bottom,  middle,  surface).  An 
analysis  was  made  of  the  freshwater  discharge  (six  canal  mouths)  versus  salinity 
at  nearby  stations.  The  expected  response  is  an  increase  in  discharge  to  be 
followed  by  a  decrease  in  measured  salinity.  Fatt  and  Wang  report  that  this  was 
not  always  the  case.  For  instance,  the  August  1985  data  showed  the  opposite 
response  for  five  of  the  six  canal  mouths.  Sudden  drops  in  salinity,  which  were 
not  preceded  by  increases  in  discharge,  are  also  reported.  These  apparently 
uncorrelated  or  counterintuitive  trends  in  salinity  values  are  of  interest.  Similar 
phenomena  were  observed  during  the  present  study.  This  is  thought  to  be  a 
consequence  of  wind-driven  circulation  or,  possibly,  localized  rainfall  and/or 
local  stratification. 
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Fatt  and  Wang  concluded  that  their  FE  salinity  transport  model  successfully 
reproduced  bay  salinity  distribution  patterns,  mixing  and  transport,  in  a  general 
way.  The  predicted  salinities  were  generally  lower  than  those  obtained  from 
field  measurements.  They  suggested  that  this  error  could  be  a  result  of  the  fact 
that  evaporation  is  not  included  in  their  model.  They  also  pointed  out  that  a 
primary  source  of  error  in  their  analysis  was  the  uncertainty  of  the  salinity  at 
model  boundaries. 

The  Wang  model  was  recently  discussed  in  the  1999  Florida  Bay  Science 
Conference,  Key  Largo,  FL.  The  model  now  includes  Barnes  Sound  and  Card 
Sound.  Calculations  have  been  made  for  salinity  and  compared  with  the  DERM 
Bay  Run  data.  However,  since  these  latest  results  have  not  yet  been  published, 
any  further  discussion  of  this  model,  including  its  underlying  assumptions  and 
efficacy  in  describing  Biscayne  Bay  hydrodynamics  and  salinity,  is  not  possible. 

Cofer-Shabica  and  Wang  (1989)  further  analyzed  the  canal  discharge  and 
salinity  data  from  the  Fatt  and  Wang  (1987)  study,  principally  that  of  the  Mowry 
Canal  (corresponding  to  the  discharge  through  S-20F  in  the  present  study).  They 
attribute  at  least  some  of  the  seemingly  unexplained  rises  and  falls  in  salinity  to 
the  discontinuous  and  highly  localized  rainfall  patterns.  Evidence  is  presented 
that  strongly  indicates  the  presence  of  fresher  water  remaining  trapped  along  the 
western  shore,  in  the  vicinity  of  the  Mowry  Canal.  Finally,  among  other  things, 
they  concluded  that  leakage  through  the  S-20F  structure  does  not  occur.  They 
suggest  it  is  reasonable  to  extend  this  conclusion  to  the  other  structures. 

A  study  by  Sengupta,  Lee,  and  Miller  (1978)  describes  a  three-dimensional 
(3-D)  finite  difference  model  of  Biscayne  Bay.  Their  model  was  developed  in 
two  versions,  one  for  sediment  transport  and  the  other  for  constituent  transport 
(without  density  coupling).  The  model  was  not  extended  to  salinity  transport. 
The  hydrodynamics  were  run  for  April  15, 1975,  and  results  compared  to  some 
statistical  norms  of  time  series  water-surface  elevation  and  surface  currents. 
Although  the  results  look  reasonable,  the  short  time  duration  and  coarse  grid 
make  it  difficult  to  learn  much  that  is  relevant  to  the  present  study. 


Features  of  present  effort 

Although  each  of  the  models  previously  discussed  has  merit,  there  are 
several  features  unique  to  the  model  presented  in  this  report  that  make  it  suitable 
for  the  evaluation  of  freshwater  impacts  on  Biscayne  Bay.  These  include: 

a.  The  inclusion  of  South  Bay  and  North  Bay  in  the  model  domain. 

b.  The  extension  of  the  ocean  boundaiy  eastward  from  the  Safety  Valve  to  a 
point  several  miles  offshore,  thereby  minimizing  the  impact  of  boundary 
effects  on  the  model  results  (especially  salinity  results). 

c.  High  grid  resolution,  especially  in  the  nearshore. 

d.  The  simultaneous  inclusion  of  groundwater  inflow,  rainfall/evaporation, 
and  wind  in  the  boundaiy  forcing. 
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e.  The  direct  application  of  freshwater  inflows  at  the  canals  mouths  (as 
opposed  to  inflows  with  salinities  calibrated  to  a  narrow  range  of  inflow 
values). 

f  The  availability  of  continuous,  yearlong,  hydrodynamic,  and  salinity  data 
at  several  locations  throughout  the  bay,  for  use  in  calibrating  and  verify¬ 
ing  the  model. 

Approach 

The  general  approach  used  in  this  numerical  effort  was  to  develop  a  compu¬ 
tational  mesh,  establish  appropriate  boundary  conditions,  and  employ  a  digital 
computer  code,  all  of  which  constitute  the  model  proper.  This  model  was  then 
used  to  compute  numerical  results  that  were  compared  to  hydrodynamic  field 
data  for  a  given  time  period  (for  this  study,  September  1997  through  October 
1997).  An  iterative  process  of  computation  of  results,  comparison  with  field 
data,  model  parameter  adjustment,  and  recalculation  and  recomparison  was  then 
undertaken.  This  cyclic  process,  referred  to  as  model  calibration,  was  continued 
until  the  model  results  and  hydrodynamic  field  data  matched  as  well  as  possible 
within  the  constraints  of  the  modeling  process.  These  constraints  refer,  pri¬ 
marily,  to  the  calibration  practice  of  adjusting  model  parameters  only  within  their 
uncertainties  or  accepted  ranges  of  variability  for  the  2-D  physics-based  numeri¬ 
cal  modeling  attempted.  Then,  the  model  was  run  for  an  additional  time  period 
(for  this  study,  November  1997  through  June  1998)  during  which  no  model 
parameters  were  adjusted.  The  results  of  this  run  were  also  compared  to  field 
data  (both  hydrodynamic  and  salinity).  This  procedure  is  referred  to  as  model 
verification.  Results  of  the  calibration  and  verification  runs  were  evaluated,  and 
additional  sensitivity  experiments  of  various  kinds  were  also  undertaken  to 
determine  the  sensitivity  of  the  model  to  uncertainties  in  various  boundary 
forcings.  Such  a  physics-based  numerical  model,  which  does  not  depend  on 
continual  adjustment  of  its  parameters,  may  then  be  run  to  assess  the  impacts  of 
various  planned  actions  with  confidence  that  the  physics  implied  by  the  planned 
action  is  responsible  for  differences  between  base  and  plan  conditions. 

The  computer  code  employed  for  this  study  is  TABS-MDS  (multidimen¬ 
sional  sediment)  (formerly  called  RMA10-WES),  a  Galerkin  based  finite  element 
algorithm  designed  to  obtain  numerical  solutions  to  one-dimensional  (1-D),  2-D, 
and/or  3-D  unsteady  free  surface  flows.  (See  King  1988;  Appendix  A  herein.) 
The  model  code  and  modeling  approach  are  out  growths  of  the  TABS-MD 
modeling  system  (Thomas  and  McAnally  1990).  The  model  boundary  conditions 
were  taken  primarily  from  data  gathered  in  Biscayne  Bay  during  the  1997  and 
1998  BNP/CHL  field  data  collection  effort  (McAdory  et  al.  2002),  Fowey  Rocks 
data  (see  http://www.noaa.gov),  and  SFWMD  canal  discharges  (McAdory  et  al. 
2002).  The  digitized  representation,  or  computational  mesh,  was  developed  early 
in  this  study  and  circulated  for  comments  among  the  study  sponsors.  Over  time, 
the  mesh  was  improved  in  response  to  comments  and  additional  data  received. 

In  conducting  this  study,  particular  emphasis  has  been  placed  on  duplicating 
the  physics  of  the  system,  rather  than  merely  trying  to  duplicate  the  system 
response,  as  given  by  the  field  data.  In  other  words,  the  range  over  which  a  given 
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model  parameter  was  adjusted  was  constrained  by  the  requirement  that  the 
parameter  value  should  represent  a  realistic  estimate  of  the  corresponding  proto¬ 
type  value  (see  prior  discussion).  This  approach  is  necessary  if  the  resulting 
model  is  to  be  a  predictive  tool.  It  minimizes  the  likelihood  of  achieving  a  good 
result  for  the  wrong  reasons,  a  potential  hazard  if  model  parameters  are  tuned  to 
force  a  desired  outcome  without  regard  to  the  physical  sensibility  of  the  param¬ 
eter  values.  Also,  focusing  on  the  physics  can  reveal  problems  with  the  initial 
modeling  assumptions  and/or  boundary  conditions  that  might  otherwise  go 
unnoticed  or  unrecognized.  As  will  be  discussed  and  argued  in  the  remainder  of 
this  report,  the  developed  2-D  numerical  model  reproduces  the  system  hydro¬ 
dynamics  and  salinity  values  and  characteristics  well,  thus  lending  itself  to  use  as 
a  physics-based,  predictive  tool  for  assessing  the  response  of  the  Biscayne  Bay 
system  to  changes  in  forcings  of  proposed  improvements  to  the  system. 
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2  Model  Description 


Grid 

Bathymetry  and  geometry 

A  previously  developed  FE  mesh  was  further  modified  for  the  purposes  of 
this  study.  The  existing  mesh  was  developed  and  circulated  among  the  study 
sponsors  early  in  the  project.  This  existing  mesh  was  generated  by  utilizing 
published  National  Oceanographic  and  Atmospheric  Administration  (NOAA) 
charts  with  bathymetry  in  feet  below  mean  lower  low  water  (mllw).  All  open 
water  regions  delineated  in  the  NOAA  charts  were  included  in  the  grid. 

Although  there  are  some  regions  along  the  western  shoreline  (particularly  in  the 
southern  part  of  the  domain)  that  have  intermittently  wetted  coastal  marshes  that 
are  not  represented  in  the  grid,  they  were  determined  to  be  either  too  small  or  too 
poorly  defined  to  be  included  in  the  grid  at  this  time.  As  these  areas  are  studied 
in  greater  detail  at  a  later  time,  these  intermittently  wetted  areas  can  be  defined 
and  included  through  either  of  two  wetting  and  drying  algorithms  available  in  the 
model  code.  Portions  of  the  current  model  already  include  the  possibility  of 
wetting  and  drying.  See  discussion  later  in  this  chapter. 

A  graphical  user  interface  (GUI),  the  Surface-Water  Modeling  System 
(SMS)  (Environmental  Modeling  Research  Laboratory  2000),  was  used  to  place 
node  points  along  the  landform  boundaries  and  interior  open  water  sections  of 
Biscayne  Bay.  SMS  is  an  elaboration  of  the  FastTABS  (BYU  1994)  visualiza¬ 
tion  system  pre-  and  post-processing  interface.  Each  node  is  identified  separately 
and  has  an  assigned  bed  elevation  at  that  location.  The  nodes  are  connected 
together  to  form  triangular  or  rectangular  elements.  The  TABS-MDS  code  uses 
an  unstructured  representation  of  the  system  geometry.  Thus,  the  placement  of 
the  nodes  need  not  be  in  a  Cartesian  arrangement.  This  unstructured  FE 
approach  then  allows  the  modeler  to  accommodate  the  irregular  shapes  of 
islands,  inlets,  and  shorelines  and  include  higher  resolution  only  where  needed. 

The  previously  developed  mesh  of  Biscayne  Bay,  used  at  the  start  of  this 
study,  contained  22,458  nodes  and  8,068  elements.  Modifications  were  made  to 
the  previously  existing  mesh.  Among  these  were:  (a)  the  inclusion  of  elements 
for  all  canals,  extending  up  to  the  structures;  (b)  finer  resolution  along  an  off¬ 
shore  reef  and  several  other  areas;  (c)  slight  relocations  and  reassignments  of 
bathymetry  to  several  nodes  along  the  western  shoreline;  (d)  the  addition  of  two 
small  channels  connecting  Manatee  Bay  and  Barnes  Sound  through  Main  Key 
and  Short  Key;  and  (e)  removal  of  a  few  minor  errors  in  node  to  node 


8 


Chapter  2  Model  Description 


connectivity.  The  resulting  FE  mesh  used  in  this  study  has  24,527  nodes  and 
8,536  elements.  For  completion,  this  mesh  is  shown  in  its  entirety  in  Figure  3. 
Figures  4,  5,  and  6  show  magnified  views  of  the  northern,  central,  and  southern 
portions  of  the  mesh,  respectively. 


Material  types 

The  TABS-MDS  model  allows  the  user  to  group  together  elements  that  share 
a  common  characteristic.  This  grouping  is  called  a  material  type.  Material  types 
provide  a  way  of  assigning  values  of  model  parameters  to  the  various  regions  of 
the  mesh  without  having  to  assign  parameters  to  each  element  individually.  The 
model  parameters  assigned  by  material  type  are  Manning’s  n  (bed  friction)  and 
the  horizontal  turbulent  mixing  parameters.  These  and  other  modeling  concepts 
are  discussed  more  in  Appendix  A.  Figure  6  shows  the  Biscayne  Bay  mesh  with 
the  major  material  types  highlighted.  Table  1  lists  these  material  types,  together 
with  a  brief  description  of  the  regions  of  the  mesh  they  represent.  Figure  7  shows 
the  material  types  of  Biscayne  Bay.  These  material  type  groupings  were  devel¬ 
oped  in  consultation  with  study  partners  and  through  use  of  the  DERM  “Bottom 
Communities  of  Biscayne  Bay”  map  (DERM  1983).'  A  description  of  the 
numerical  values  of  parameters  assigned  to  the  model  by  material  type  is 
presented  later  in  this  chapter. 


Table  1 

Major  Material  Types  and  Descriptions 

Materia!  Type 

Description 

1 

Offshore,  smooth  bed 

2 

Middle  Bay,  grass  beds 

3 

Card  Bank,  grass  beds 

4 

South  Bay,  inlets 

5 

Offshore  reefs 

6 

Little  Card  Sound,  grass  beds  _J 

7 

Card  Sound,  grass  beds 

8 

North  Bay,  inlets  (Including  Government  Cut  and  Bakers  Haulover  Inlet) 

10 

Canals 

11 

Middle  Bay,  western  shoreline 

12 

Barnes  Sound  and  Manatee  Bay,  grass  beds 

13 

Cutter  Bank 

14 

South  and  Southeast  Bay,  southwest  of  Elliott  Key 

15 

Featherbed  Bank,  Ragged  Keys,  and  Soldier  Key 

16 

North  Bay,  west  of  Virginia  Key 

17 

North  Bay,  Julia  Tuttle  Causeway  to  Bakers  Haulover  Inlet 

18 

North  Bay,  South  of  Julia  Tuttle  Causeway 

19 

Shoals,  Julia  Tuttle  Causeway  to  Bakers  Haulover  Inlet 

20 

Offshore  Boundary 

1  Miami-Dade  Department  of  Environmental  Resources  Management.  (1983).  “Bottom 
communities  of  Biscayne  Bay”  Map,  Published  in  Cooperation  with:  State  of  Florida  Department 
of  Environmental  Regulation,  United  States  Department  of  the  Interior,  National  Park  Service, 
Miami,  FL. 
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Resolution 


The  surface  areas  of  the  elements  in  the  mesh  vary  over  several  orders  of 
magnitude,  from  as  small  as  1  x  104  ft2  to  as  large  as  3  x  107  ft2.  Locations  where 
velocity  and/or  salinity  gradients  are  expected  to  be  high  are  highly  resolved,  as 
are  regions  of  special  interest  to  the  study  (such  as  the  western  shore).  In  regions 
where  gradients  are  expected  to  be  small,  or  where  there  is  less  interest  in  resolv¬ 
ing  the  local  current/salinity  field  (such  as  in  the  offshore),  the  resolution  is 
coarser.  The  unstructured  nature  of  the  finite  element  mesh  allows  this  approach 
to  be  taken.  The  mesh  for  Biscayne  Bay  was  developed  by  taking  full  advantage 
of  this  capability. 


Model  Boundary  Conditions  and  Initial  Conditions 

The  TABS-MDS  model  requires  initial  and  boundary  conditions,  as  well  as  a 
number  of  modeling  parameters,  to  be  specified  via  input  files.  The  governing 
equations  for  hydrodynamics  require  that  either  water-surface  elevation  or 
velocity  (or  flow  rate)  be  specified  at  all  nodes  along  open  boundaries  (i.e., 
boundaries  other  than  land  boundaries).  In  addition  to  these  necessary  condi¬ 
tions,  the  Biscayne  Bay  model  was  driven  with  additional  inputs,  including 
applied  wind  stress,  rainfall/evaporation,  and  applied  groundwater.  At  the  land 
boundaries,  TABS-MDS  imposes  a  slip  condition,  forcing  the  velocity  to  be 
parallel  to  the  shoreline.  That  is,  the  velocity  magnitude  is  solved  for,  while  the 
velocity  direction  is  forced  to  be  tangent  to  the  land  boundary. 

To  solve  for  salt  transport,  it  is  necessary  to  impose  salinity  boundary  condi¬ 
tions.  For  Biscayne  Bay,  a  salt  concentration  was  specified  along  the  offshore 
boundary  for  all  flow  into  the  model  domain.  The  concentration  of  flow  exiting 
the  model  domain  through  the  offshore  boundary  was  not  specified.  (To  specify 
the  concentration  of  flow  exiting  the  model  in  this  way  would  be  physically 
nonsensical  because  the  salinity  of  the  near  boundary  water  that  will  be  pushed 
out  of  the  ocean  boundary  is  intended  to  be  determined  by  the  internal  dynamics 
of  the  model.)  The  concentration  was  specified  for  all  inflows  to  the  system,  as 
well  as  for  rainfall/evaporation  and  groundwater. 

The  unsteady  nature  of  the  Biscayne  Bay  problem  necessitates  that  a  time¬ 
stepping  procedure  be  used,  so  that  the  boundary  conditions  can  be  changed  for 
each  new  time-step.  Hence,  for  Biscayne  Bay,  all  the  applied  boundary  condi¬ 
tions  were  given  as  time  series.  In  this  way  a  simulation  of  the  entire  bay  was 
performed.  The  influences  of  changes  in  the  applied  boundary  conditions  with 
respect  to  time  were  propagated  throughout  the  model  domain  via  the  numerical 
solution  of  the  governing  equations.  Removing  these  boundaries,  particularly  the 
ocean  salinity  and  tidal  boundaries,  far  from  the  area  of  interest  (i.e.,  the  Bay 
proper)  minimizes  the  pegging  of  bay  salinities  and  water  levels  by  the  specified 
boundary  values.  This  allows  the  values  of  salinity  and  water  level,  for  instance, 
in  the  area  of  interest  to  be  shaped  by  the  modeled  physics.  Initial  conditions  are 
also  necessary  for  the  modeled  physical  parameters.  To  assign  initial  conditions 
for  the  model,  an  initial  solution  field  is  generated  for  both  the  hydrodynamics 
and  the  salinity. 
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The  entire  set  of  boundary  conditions  assembled  for  the  calculations  dis¬ 
cussed  in  this  report  was  12  months  long,  spanning  the  time  from  August  1997 
through  July  1998.  However,  no  rainfall  or  evaporation  data  were  available  for 
August  1997  or  July  1998  (rainfall  data  were  also  unavailable  for  September 
1997,  but  for  this  month  an  estimate  of  the  rainfall  was  taken  from  raw  pan 
evaporation  data,  which  record  net  precipitation).  Since  rainfall  and  evaporation 
are  considered  significant  parameters  in  the  Biscayne  Bay  system,  it  was  not 
considered  suitable  to  include  these  months  in  the  simulation  without  these  data. 
Therefore,  the  entire  period  of  simulation  consists  of  the  10-month  period  extend¬ 
ing  from  September  1997  through  June  1998.  The  hydrodynamic  model  was 
calibrated  with  the  September  through  October  1997  data,  and  verified  against 
the  November  1997  through  June  1998  data.  No  calibration  was  figured  for  the 
salinity  (see  Chapter  3).  The  salinity  was  permitted  to  spin  up  through 
September  1997  and  was  verified  against  the  October  1997  through  June  1998 
data. 


Offshore  tidal  boundary 

The  time  series  of  water-surface  elevations  used  along  the  offshore  boundary 
was  extrapolated  from  measured  water-surface  elevations  at  Gauge  9  conductiv¬ 
ity,  temperature,  and  depth  recorder  No.  9  (CTD9)  of  the  BNP/CHL  data  set. 

The  signal  from  CTD9  was  modified  such  that  it  matched  the  measured  signal  at 
the  Virginia  Key  National  Ocean  Service  (NOS)  station  (No.  8723214),  located 
at  the  Rosenstiel  School  of  Marine  and  Atmospheric  Science  (RSMAS)  for  the 
calibration  time  interval.  The  Virginia  Key  station  was  used  because  the 
National  Geodetics  vertical  datum  (NGVD29)  is  well  established  at  this  station. 
The  Virginia  Key  station  data  were  not  used  in  place  of  the  CTD9  water-surface 
elevation  data  because  the  CTD9  data  provided  the  more  complete  data  set  for 
the  entire  data  collection  period.  The  following  adjustments  were  made  to  the 
CDT9  data  to  establish  a  match  with  the  Virginia  Key  data: 

Phase  shift  =  -  45  min 

Vertical  offset  =  +0.7  ft  (including  0.67  ft  for  adjustment  from  mllw  to 
NGVD29) 

Signal  amplification  =  27  percent  (amplitude  multiplied  by  1.27) 

The  phase  shift  accounts  for  the  time  of  travel  from  the  model  boundary  to 
the  modeled  position  of  the  Virginia  Key  station.  The  vertical  offset  ensures  the 
model  data  are  appropriate  at  the  Virginia  Key  station.  And  the  signal  amplifica¬ 
tion  accounts  for  deamplification  of  the  tidal  signal  as  it  moves  toward  the  coast. 
It  was  found  that  by  making  these  adjustments  to  the  CTD9  data,  TABS-MDS 
consistently  reproduced  the  tide  signal  at  the  NOS  Virginia  Key  station  through 
the  calibration  period.  The  model  and  field  data  for  the  calibration  period  are 
shown  for  water-surface  elevation  in  Plates  1-15.  Virginia  Key  data  are  given  in 
Plate  15.  The  time  series  generated  for  the  model  is  given  in  30-min  increments, 
which  corresponds  to  the  length  of  the  time-step  used  in  the  simulation. 
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Offshore  salinity 


The  salinity  along  the  ocean  boundary  was  specified  as  a  constant  value  for 
each  month  within  the  data  set.  To  establish  the  monthly  values  to  be  applied  to 
the  offshore  boundary,  data  were  examined  from  two  different  locations:  the 
Fowey  Rocks  Lighthouse,  located  about  414  miles  east  of  Soldier  Key,  and  the 
BNP  gauge  at  Alina's  Reef,  located  about  3  miles  east  of  Card  Sound  (Figure  8). 
Both  data  sets  are  shown  in  Figure  9.  Examination  of  the  data  reveals  trends  that 
appear  to  indicate  instrument  fouling  at  the  Fowey  Rocks  gauge.  The  evidence 
for  this  is  the  periodic  appearance  of  a  gradual  downward  trend  in  salinity, 
followed  by  a  sudden  increase  that  corresponds,  in  time,  to  the  maintenance  of 
the  gauge.  Such  behavior  is  indicative  of  gauge  fouling.  Because  of  this 
apparent  fouling  of  the  Fowey  Rocks  data,  the  monthly  averaged  Alina’s  Reef 
data  were  used  as  the  ocean  salinity  boundary  condition  (Figure  9  and  Table  2). 


Wind 

The  wind  applied  to  the  model  was  based 
on  data  from  the  BNP/CHL  gauge  located  at 
the  Convoy  Point  station.  Wind  values  were 
applied  over  the  entire  bay  portion  of  the 
computational  domain,  with  the  value  updated 
at  30-min.  intervals.  Wind  values  were  not 
applied  over  the  ocean  portions  of  the  model 
domain.  Qualitative  comparisons  of  the 
Convoy  Point  wind  data  with  wind  data  taken 
at  the  NOAA  gauge  at  Fowey  Rocks  suggest 
that  the  dominant  wind  direction  is  relatively 
homogeneous  in  the  area  of  the  bay,  with  the 
biggest  overall  difference  being  a  decrease  in 
velocity  magnitude  from  the  Fowey  Rocks 
station  to  the  Convoy  Point  station.  Smith 
also  reached  this  conclusion  and  provided 
quantitative  comparisons,  finding  the  Fowey  Rocks-Convoy  Point  velocity 
magnitudes  to  be  in  the  ratio  1.45  (McAdory  et  al.  2002).  The  Convoy  Point 
wind  speed  data  for  the  entire  simulation  period  are  shown  in  Figure  10. 

The  raw  Convoy  Point  data  were  manipulated  to  account  for  differences  in 
magnitude  due  to  height  of  collection  and  the  near  land  location  of  the  Convoy 
Point  station.  The  raw  data  were  corrected  to  the  standard  10-m  height  above  the 
water  using  the  factor  1.16  (Hsu  1 988).  These  height  corrected  data  were  then 
used  to  calculate  wind  stress  on  the  water  surface  in  North  Bay,  South  Bay  (Card 
and  Bames  Sounds),  and  along  the  western  shoreline  in  Central  Bay  using  the 
formulation  of  Wu  (1980).  This  application  was  made  assuming  these  areas  are 
relatively  sheltered  as  compared  with  the  main  Central  Bay  area.  To  account  for 
higher  wind  velocity  magnitudes  over  the  open  water  of  Central  Bay,  half  the 
factor  relating  the  Fowey  Rocks  ocean  wind  data  to  the  Convoy  Point  nearshore 
wind  data  was  applied,  namely  1.225  (Figure  1 1). 


Table  2 

Salinity  at  Fowey  Rocks  and 
Alina’s  Reef 

Month 

Alina’s  Reef 
Average  (all 
data) 

Fowey  Rocks 
Average  (all 
data) 

Sept.  97 

34.9 

35.1 

Oct.  97 

35.4 

35.8 

Nov.  97 

35.7 

35.2 

Dec.  97 

35.4 

35.0 

Jan.  98 

35.9 

34.6 

Feb.  98 

35.7 

33.7 

Mar.  98 

35.9 

33.8 

Apr  98 

36.6 

34.7 

May  98 

36.6 

35.8 

June  98 

34.9 

35.2 
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Rainfall  and  evaporation 

Rainfall  and  evaporation  data  were  taken  from  the  weather  station  at  Convoy 
Point.  These  values  were  applied  over  the  entire  model  domain,  with  the  values 
updated  at  30-min.  intervals.  The  spatial  distribution  of  rainfall  and  evaporation 
in  Biscayne  Bay  can  be  significantly  heterogeneous,  especially  in  the  summer 
months,  when  storms  tend  to  be  intense  and  localized.  Use  of  a  set  of  uniform 
rainfall  and  evaporation  data  from  a  station  located  in  the  system,  e.g.,  the 
Convoy  Point  station,  means  that  the  data  are  assumed,  in  the  large,  to  be  repre¬ 
sentative  of  the  system.  Spatially  varying  data  can  be  incorporated  into  the 
model,  when  available,  as  desired.  Figures  12  and  13  show  rainfall  and  evapo¬ 
ration  for  the  entire  simulation  period. 


Inflows:  Canals  and  rivers 

Freshwater  inflows  from  the  various  canals  and  rivers  that  flow  into  Biscayne 
Bay  are  regulated  by  16  structures.  Figure  8  shows  the  locations  of  these  struc¬ 
tures.  These  structures  are  managed  by  the  SFWMD  (Van  Horn  1996).  Rating 
curves,  developed  by  SFWMD  and  the  U.  S.  Army  Corps  of  Engineers 
(US  ACE),  exist  for  each  of  these  structures  (Swain  et  al.  1997),  and  revised 
curves  were  developed  by  Swain  et  al.  (1997).  Since  the  integration  of  these  new 
curves  into  the  SFWMD  canal  flow  data  sets  was  not  complete  at  the  time  of  this 
study,  sponsors  recommended  use  of  the  earlier  rating  curves  (Imru  2001).  In 
any  case,  flow  data  used  in  the  modeling  effort  were  flows  provided  to  CHL  by 
the  SFWMD  (McAdory  et  al.  2002).  These  flows  were  generated  by  SFWMD 
from  the  DBHYDRO  database  (available  from  SFWMD)  using  the  existing 
rating  curves  to  generated  flow  hydrographs  for  each  structure.  For  use  in  the 
modeling,  the  hydrographs  were  averaged  over  30-min.  intervals,  and  the  result¬ 
ing  values  were  applied  as  inflows  to  the  model.  Checks  were  made  to  ensure 
that  this  30-min  averaging  did  not  change  the  net  flow  into  the  bay. 

There  is  significant  uncertainty  concerning  the  accuracy  of  the  rating  curves, 
particularly  for  extreme  events.  The  coefficients  used  for  these  curves  are  based 
on  logarithmic  regression  analyses  of  data  that  represent  limited  ranges  of  flow 
(Swain  et  al.  1997).  Since  several  of  the  events  that  occur  during  the  period  of 
record  exceed  these  ranges,  it  is  uncertain  whether  or  not  the  coefficients  can  be 
trusted  in  these  high  flow  cases.  This  question  is  the  subject  of  ongoing  investi¬ 
gation  at  the  SFWMD.  In  a  recent  report,  personnel  at  the  SFWMD  supply  esti¬ 
mates  of  the  absolute  percent  error  in  the  prediction  of  the  flow  at  each  of  the 
structures  (Imru  2001).  The  errors  reported  for  the  SFWMD/USACE  rating 
curves  are  given  in  Table  3. 

The  hydrographs  of  each  of  the  16  structures  for  the  entire  simulation  period, 
as  used  in  the  modeling  effort,  are  given  in  Plates  1  through  8.  The  salinity  of 
these  inflows  was  specified  as  0.0.  Although  some  small  values  of  salinity  may 
occur  upstream  of  the  structures,  it  was  assumed  that  flow  of  any  significant 
duration  would  flush  this  trace  salinity  from  the  canals.  Finally,  negative  flow 
values  were  set  to  zero  in  the  modeling  effort.  Since  the  canal  structures  are 
specifically  designed  to  avoid  such  flows,  these  apparent  flows  are  considered  to 
be  an  artifact  of  the  telemetric  process. 
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Table  3 

Estimated  Error  in  Flow  Prediction  at  Biscayne  Bay  Structures 

Salinity  Control  Structure 

Absolute  Percent  Error  In  Flow  Estimate  for  II 
SFWMD/COE  Rating  Curves 

S-29 

19 

G-58 

No  data 

S-28 

29 

S-27 

23 

S-26 

38 

S-25B 

9 

S-25 

No  data 

G-93 

20 

S-22 

9 

S-123 

12 

S-21 

24 

S-21A 

22 

S-20G 

29 

S-20F 

50 

S-20 

24 

S-197 

17 

Inflows:  Jewfish  Creek 

Jewfish  Creek  is  a  small  inlet  connecting  Bames  Sound  to  Blackwater 
Sound.  For  this  model  study,  it  was  treated  as  a  flow  boundary.  Flow  and 
salinity  data  were  available  for  Jewfish  Creek  throughout  the  period  of  simu¬ 
lation.  However,  intermittent  gaps  in  the  data  record  exist.  To  fill  these  gaps, 
synthetic  data  were  generated,  based  on  a  simple  correlation  between  water- 
surface  elevation  data  from  Manatee  Bay  (CTD1)  and  the  flow  at  Jewfish  Creek 
for  periods  when  flow  data  exist.  Since  the  volume  of  flow  passing  through 
Jewfish  Creek  is  small  relative  to  the  tidal  exchange  in  Bames  Sound,  and  since 
the  salinity  is  similar  to  the  salinity  found  in  Bames  Sound,  the  error  introduced 
into  the  hydrodynamic  and  salinity  modeling  by  using  synthetic  data  to  fill  gaps 
in  the  data  record  is  assumed  to  be  negligible.  For  reasons  similar  to  those 
discussed  for  the  offshore  salinity  boundary,  salinity  for  inflow,  only,  into  Bames 
Sound  through  Jewfish  Creek  was  specified.  These  salinity  values  were  based  on 
field  data  taken  at  Jewfish  Creek.  Figure  14  shows  the  Jewfish  Creek  hydrograph 
and  salinity  for  the  entire  simulation  period  as  used  in  the  modeling  effort. 


Inflows:  Groundwater 

The  groundwater  inflow  was  taken  directly  from  model  results  generated  by 
the  U.  S.  Geological  Survey  (USGS)  (Langevin  2001).  The  groundwater  inflow 
was  applied  in  the  following  locations:  along  the  western  shoreline  of  Biscayne 
Bay  north  of  structure  S-123  (this  includes  approximately  21  miles  of  shoreline); 
along  the  western  shoreline  of  Biscayne  Bay  south  of  structure  S-123  (this 
includes  approximately  26  miles  of  shoreline);  in  the  Miami  River;  and  in  the 
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Coral  Gables  Canal.  Figure  15  shows  these  locations.  The  groundwater  inflow 
was  distributed  as  follows: 

a.  37.5  percent  was  applied  along  the  shoreline  north  of  S-123. 

b.  12.5  percent  was  applied  along  the  shoreline  south  of  S-123. 

c.  25.0  percent  was  applied  in  the  Miami  River. 

d.  25  .0  percent  was  applied  in  the  Coral  Gables  Canal. 

No  groundwater  was  applied  during  September  1997.  However,  September 
1997  is  a  spin-up  month  for  the  model  salinity;  this  omission  was  assumed  to 
have  negligible  impact  on  the  modeling  results  for  subsequent  months  of  simu¬ 
lating.  This  assumption  is  reinforced  by  the  small  amount  of  groundwater  inflow 
when  compared  with  the  September  1997  canal  inflows.  The  groundwater  was 
assigned  a  salinity  of  17.5,  again  based  on  the  model  results  of  the  USGS 
(Langevin  2001).  The  groundwater  inflows  applied  to  the  model  for  the  entire 
simulation  period  (except  September  1997)  are  given  in  Figure  16. 


Offshore  currents 

Offshore  currents  are  considered  an  important  mechanism  for  sweeping  away 
low  salinity  water  from  Biscayne  Bay.  To  simulate  these  offshore  sweeping 
currents,  and  thus  allow  for  their  effects  to  be  realized  in  the  area  of  interest  (i.e., 
in  the  bay),  a  north-south  slope  was  placed  on  the  ocean  boundary  of  the  compu¬ 
tational  domain.  This  slope  generates  offshore  currents  that  simulate  the  effects, 
in  the  bay,  of  the  actual  offshore  currents.  This  is  particularly  the  case  in  terms 
of  sweeping  Biscayne  Bay  water  that  exits  the  bay  away  from  its  exit  point.  The 
sweeping-current-inducing  slope  of  the  offshore  area  was  varied  in  time  accord¬ 
ing  to  the  direction  of  the  winds.  The  induced  sweeping  currents  thus  vary  from 
northerly  to  southerly  in  direction,  depending  on  the  offshore  wind  direction. 

The  predominant  direction,  as  evidenced  by  the  wind,  is  toward  the  north.  The 
result  is  a  sweeping  current  with  the  proper  order  of  magnitude  for  both  current 
direction  and  current  magnitude.  This  result  is  sufficient  to  account,  in  the  bay, 
for  sweeping  effects  outside  of  the  bay  at  the  level  of  this  modeling  effort. 


Initial  conditions 

The  initial  conditions  for  the  hydrodynamics  were  generated  by  the  model. 
The  model  was  used  to  solve  for  a  steady  state  solution  for  Biscayne  Bay,  and 
this  solution  was  used  as  the  initial  condition  field  for  velocities  and  depth.  The 
initial  conditions  for  salinity  were  interpolated  from  field  data.  Salinity  values 
were  extracted  from  the  BNP/CHL  data  set,  and  these  were  used,  together  with  a 
2-D  interpolation  scheme,  to  generate  an  initial  salinity  field  for  the  entire  model 
domain. 
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Computer  Code 

General  discussion 

The  numerical  model  used  for  this  study  is  TABS-MDS.  TABS-MDS  is  a 
Galerkin-based  FE  solution  for  unsteady  open  channel  flow  and  sediment 
transport.  The  model  is  based  on  RMA10,  a  code  originally  developed  by  Ian 
King  of  Resource  Management  Associates  (King  1988)  under  contract  to  the 
Corps.  The  CHL  version  of  the  model  has  been  modified  to  such  an  extent  that, 
to  avoid  confusion,  the  CHL  version  has  been  assigned  the  name  TABS-MDS. 

The  code  is  capable  of  simulating  hydrodynamics  in  1-D,  2-D  and/or  3-D 
modes.  It  solves  conservation  of  fluid  mass  equations,  horizontal  momentum 
equations,  and  equations  of  state  (for  the  determination  of  density  as  a  function  of 
salinity,  temperature,  and/or  suspended  sediment),  and  convection/diffusion 
equations.  To  improve  computational  efficiency,  the  hydrostatic  assumption  is 
implemented.  This  assumption  implies  that  vertical  accelerations  are  negligible. 
This  approximation  is  valid,  and  widely  used,  in  estuarine  application  of  model 
codes. 

TABS-MDS  is  both  unstructured  and  implicit.  An  unstructured  code  maxi¬ 
mizes  the  ability  to  resolve  complicated  shoreline  or  other  geometry,  and  an 
implicit  code  allows  the  time-step  length  to  be  largely  unconstrained  by  the  grid 
resolution  and  local  velocities.  These  are  generally  needed  for  generating 
accurate  representations  of  estuarine  systems  with  complicated  shorelines  and 
dramatic  spatial  differences  in  velocity,  such  as  is  found  in  Biscayne  Bay.  For 
more  details  concerning  the  TABS-MDS  model,  consult  Appendix  A. 


Density  coupling 

For  this  study,  TABS-MDS  was  run  in  2-D  vertically  averaged  mode.  By 
definition,  this  eliminates  the  possibility  of  modeling  salinity  or  velocity  strati¬ 
fied  profiles  in  the  bay.  However,  coupling  the  density  in  a  2-D  calculation  does 
permit  the  inclusion  of  horizontal  density  gradient  effects,  which  are  the  domi¬ 
nant  density  gradient  contribution  in  well-mixed  systems.  This  density  coupling 
was  enabled  for  this  study. 


Modeling  of  intermittently  wetted  regions 

TABS-MDS  is  equipped  with  a  feature  that  permits  the  water  and  salt  mass 
storage  in  intermittently  wetted  regions  of  the  domain  (e.g.,  marshes,  shoals,  etc.) 
to  be  accounted  for  in  the  calculation.  This  feature  is  called  marsh  porosity,  as 
was  originally  proposed  by  Roig  (1995).  The  marsh  porosity  method  distributes 
the  water  in  a  partially  submerged  element  over  the  entire  element,  thereby 
creating  a  thin  layer  of  water  submerging  100  percent  of  the  surface  area  of  the 
element.  This  method  conserves  both  mass  and  momentum.  Since  the  distribu¬ 
tion  of  water  over  each  element  is  likely  different  than  the  distribution  observed 
in  the  field,  the  friction  loss  across  a  marsh  porosity  element  may  be  different  in 
the  model  than  in  the  field.  Therefore,  if  marsh  porosity  elements  are  located  in  a 
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region  of  the  mesh  where  momentum  conservation  is  important,  care  must  be 
exercised  in  the  selection  and  calibration  of  friction  parameters  for  these 
elements.  For  more  information  on  marsh  porosity,  consult  Roig  (1995). 

For  Biscayne  Bay,  marsh  porosity  was  activated  in  the  mesh  at  any  node  for 
which  the  depth  is  0.5  ft  or  less.  There  are  several  areas  of  the  mesh  where  this  is 
a  possibility.  The  shallowest  depths  in  Biscayne  Bay  occur  at  Card  Bank  and  the 
shoals  in  the  North  Bay  (material  types  3  and  19  (Table  1  and  Figure  7)).  The 
bed  elevations  in  these  regions  get  as  high  as  -0.8  ft,  referenced  to  mllw.  There¬ 
fore,  any  time  the  water-surface  elevation  (relative  to  mllw)  is  less  than  -9  cm 
(-0.3  ft)  in  these  regions,  marsh  porosity  will  become  active  (i.e.,  these  regions 
will  be  effectively  dry).  This  means  that,  although  marsh  porosity  may  be  inter¬ 
mittently  activated  in  these  regions,  it  is  never  active  for  an  extended  duration. 

As  mentioned  earlier  in  this  chapter,  further  application  of  marsh  porosity 
wetting  and  drying  can  be  added  to  other  areas  of  the  modeled  bay,  such  as  the 
western  areas,  as  desired  for  future  model  applications. 


Bottom  friction 

Bottom  friction  in  TABS-MDS  is  represented  with  a  modified  form  of 
Manning’s  equation.  This  modification  is  intended  to  correct  for  depth 
dependence  in  the  calculation  of  bottom  friction.  For  more  on  this  method, 
consult  Appendix  A. 


Horizontal  turbulent  mixing  and  diffusion 

Horizontal  turbulent  mixing  and  diffusion  are  represented  in  TABS-MDS  by 
the  method  of  Smagorinsky  (1963).  This  method,  originally  developed  for 
atmospheric  models,  has  been  widely  used  in  hydrodynamic  models  (Hamrick 
1992).  Given  the  appropriate  choice  of  coefficients,  this  method  has  demon¬ 
strated  the  capacity  to  model  sub-grid  scale  dissipation  in  hydrodynamic  systems 
(Speziale  1998).  Details  of  the  implementation  of  the  Smagorinsky  method  into 
the  Biscayne  Bay  model  are  given  in  the  calibration  section  of  Chapter  3.  For 
more  details  on  the  Smagorinsky  method,  consult  Appendix  A. 


Wind  stress 

The  Wu  formulation  (Wu  1980)  is  used  to  calculate  wind  stress  on  the  bay 
water  surface.  This  method  is  assumed  to  be  adequate  given  the  estimates  of 
wind  velocity  used. 
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3  Model  Calibration  and 
Verification 


Overview 

Calibration  refers  to  the  process  of  running  a  simulation  for  a  given  period  of 
record  and  adjusting  selected  parameters  to  match  the  model  results  to  the  field 
data.  Verification  refers  to  the  process  of  running  a  simulation  for  a  period  of 
record  other  than  the  one  used  for  calibration  and  comparing  model  simulation 
results  and  field  data  without  adjusting  any  parameters.  As  previously  stated  in 
Chapter  1,  the  hydrodynamics  were  calibrated  with  data  from  September  through 
October  1997,  and  verified  against  data  from  November  1997  through  June  1998. 
The  salinity  was  not  calibrated;  that  is,  no  further  adjustments  were  made  to 
model  parameters,  beyond  what  adjustments  were  made  for  the  hydrodynamic 
calibration,  to  alter  the  salinity  results.  Salinity  calculations  were  permitted  to 
spin  up  through  September  1997,  and  then  the  salinity  was  verified  against  data 
from  October  1997  through  June  1998. 


Prototype  Data 

The  following  is  a  brief  discussion  of  the  prototype  data.  For  more  detailed 
information,  consult  McAdory  et  al.  (2002). 


Biscayne  National  Park/Coastal  and  Hydraulics  Laboratory  data 

This  is  an  extensive  data  set,  and  the  main  one  used  for  this  study.  It  consists 
of  four  major  parts:  water-surface  elevation  data,  salinity  data,  velocity  data,  and 
flow,  or  discharge,  data. 

Water-surface  elevation  data.  These  data  were  collected  at  12  CTD  moni¬ 
toring  gauges  located  throughout  the  model  domain.  The  gauges  were  attached 
to  the  bed  and  floated  in  the  water  column  with  submerged  buoys,  which  placed 
the  sensors  approximately  1  ft  from  the  bed.  Data  were  collected  synoptically  at 
1 1  CTD  locations  from  June  1997  to  June  1998  and  from  January  to  June  1998 
for  the  CTD  10  in  the  mouth  of  the  Miami  River.  A  data  record  was  recorded 
every  15  min.  Each  data  record  consisted  of  time,  temperature,  depth,  pH, 
dissolved  oxygen,  and  conductivity.  Figure  17  shows  the  locations  of  these 
gauges. 
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Salinity  data.  Salinity  data  were  calculated  internally  to  each  CTD  from  the 
measured  conductivity  and  temperature  values  at  each  CTD  and  stored  internally 
for  later  retrieval  for  the  same  period  of  record  that  the  water-surface  elevation 
data  were  collected.  Therefore,  salinity  data  were  also  collected  approximately 
1  ft  from  the  bed.  Water-surface  elevation  and  salinity  values  were  recorded 
independently,  however,  so  data  gaps  in  one  of  these  quantities  may  not  be 
present  in  the  other. 

Velocity  data.  Velocity  data  were  collected  at  five  locations  in  the  bay, 
using  bottom  mounted  Acoustic  Doppler  Profiling  (ADP)  instruments.  At  each 
gauge,  the  velocities  were  reported  at  multiple  depths.  This  was  done  by  divid¬ 
ing  the  water  column  into  1 1.8  in.  sections,  or  “bins,”  and  measuring  the  velocity 
in  each  bin  (note  that  the  bottom  23.6  in.  of  the  water  column  are  unmeasured 
because  of  limitations  of  the  instrument).  Measurements  were  taken  every 
15  min.  Figure  17  shows  the  locations  of  these  instruments. 

Note  that,  although  velocity  data  were  collected  at  ADP4,  the  depth  sensor 
was  not  recording  properly,  and  therefore  the  locations  of  each  velocity  measure¬ 
ment  within  the  water  column  were  impossible  to  determine.  Since  the  method 
of  depth  averaging  employed  for  this  study  depends  on  knowing  the  locations 
within  the  water  column  of  each  velocity  measurement,  it  was  decided  to  omit 
the  data  from  ADP4  from  consideration  in  this  study. 

Flow  data.  Flow  data  were  collected  with  Acoustic  Doppler  Current 
Profilers  (ADCP)  mounted  aboard  vessels.  Two  separate  data  sets  were 
gathered,  one  in  October  and  one  in  February.  Each  data  set  consists  of  27  tran¬ 
sects,  each  sampled  for  approximately  9  hr.  Figurel8  depicts  the  locations  of 
these  transects.  Note  that  Transect  24  was  not  taken  in  the  February  collection 
period. 


DERM  “Bay  Run”  data 

The  DERM  Bay  Run  water  quality  data  are  salinity  data  (and  a  number  of 
other  parameters)  collected  over  multiple  sites  in  Biscayne  Bay  (DERM  various). 
One  measurement  is  taken  per  month,  and  each  measurement  consists  of  at  least 
two  readings  at  different  depths  in  the  water  column,  yielding  an  instantaneous 
local  salinity  profile.  The  readings  always  consist  of  at  least  one  reading  at  the 
water  surface,  and  one  reading  at  a  depth  of  3.3  ft.  If  the  local  water  depth  at  the 
sampling  site  is  greater  than  3.6  ft  at  the  time  of  the  reading,  then  a  third  sample 
is  taken  at  the  bottom.  The  locations  of  these  sampling  points  are  given  in 
Figures  19  through  21. 


Sea-bird  data 

Periodically,  the  BNP/CHL  gauges  were  serviced  and  replaced.  Each  time 
this  occurred  for  the  CTD  locations,  a  hand-held  salinity  gauge  was  lowered  over 
the  side  of  the  vessel,  and  a  salinity  profile  was  recorded.  Because  of  limitations 
of  the  instrument,  however,  this  profile  never  includes  a  salinity  measurement  in 
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the  upper  3.3  ft  of  the  water  column.  A  Sea-Bird  instrument  was  used  for  this 
profiling  (McAdoiy  2002). 


Model  Calibration 

The  hydrodynamic  model  was  calibrated  over  the  time  period  from 
September  through  October  1997.  The  model  was  run  with  a  time-step  of  0.5  hr, 
giving  approximately  25  time-steps  per  tidal  cycle.  The  tide  signals  at  the  12 
CTD  locations  were  compared  to  the  model  results  at  the  CTD  locations  in  the 
model,  and  adjustments  were  made  to  the  friction  parameters  in  the  model. 
Velocity  and  discharge  calculations  were  then  compared  to  field  data  and  further, 
more  refined  adjustments  were  made  to  the  model  friction. 

The  eddy  viscosities  in  the  model  were  adjusted  during  the  period  of  calibra¬ 
tion,  but  only  for  the  purposes  of  achieving  numerical  stability.  As  was  stated  in 
Chapter  2,  the  eddy  viscosity  in  the  Biscayne  Bay  study  was  controlled  by  the 
Smagorinsky  method  (Smagorinsky  1963).  In  TABS-MDS,  there  are  actually 
two  quantities  that  control  the  eddy  viscosity  when  using  the  Smagorinsky 
method.  The  first  quantity  is  the  Smagorinsky  parameter,  which  is  adjusted  with 
the  Smagorinsky  coefficient.  The  second  is  a  minimum  allowable  eddy  viscosity 
parameter.  This  value  is  used  to  supply  the  model  with  sufficient  diffusion  to 
ensure  stability.  At  each  computational  point  and  at  each  time-step,  the  model 
calculates  both  values  of  eddy  viscosity  (one  with  the  Smagorinsky  parameter 
and  the  other  with  the  minimum  allowable  eddy  viscosity  parameter)  and  selects 
the  maximum  of  the  two  for  the  value  it  uses.  A  value  of  0.05  was  selected  for 
the  Smagorinsky  coefficient,  which  is  in  the  range  of  accepted  values  typically 
used  in  estuarine  systems  (Speziale  1998).  Then,  a  nominal  initial  estimate  of  the 
minimum  eddy  viscosity  parameter  was  selected.  The  model  was  allowed  to  run, 
and  local  instabilities  were  eliminated  by  adjusting  the  minimum  eddy  viscosity 
parameter  until  the  instabilities  were  corrected.  Since  the  minimum  eddy  vis¬ 
cosity  parameter  can  be  applied  by  material  type,  this  adjustment  could  be  done 
locally,  without  elevating  the  global  value  of  the  minimum  eddy  viscosity.  In 
some  cases,  when  only  a  few  elements  were  manifesting  the  instability,  a  new 
material  type  was  created  to  allow  the  minimum  eddy  viscosity  parameter  to  be 
adjusted  for  those  elements  only.  This  method  of  adjustment  ensured  that  the 
minimum  eddy  viscosity  that  was  applied  to  the  model  was  just  sufficient  to 
prevent  numerical  instabilities  from  occurring. 

Although  the  wind  speed/stress  was  not  calibrated,  the  model  and  field  tide 
and  velocity  data  for  October  were  filtered  to  remove  the  tidal  signal,  and  the 
remaining  wind-induced  signals  were  compared.  This  is  discussed  further  in  the 
“Qualitative  Analysis”  section  later  in  this  chapter. 

No  calibration  was  performed  for  the  salinity  modeling.  The  horizontal 
diffusion  coefficients  were  determined  by  the  same  process  used  to  adjust  the 
eddy  viscosity.  That  is,  a  value  of  0.05  was  selected  for  the  Smagorinsky  coeffi¬ 
cient,  a  nominal  minimum  diffusion  parameter  was  selected,  and  the  minimum 
diffusion  parameter  was  adjusted  until  any  numerical  instabilities  were  elimi¬ 
nated.  Hence,  as  with  the  eddy  viscosity,  the  only  adjustments  applied  to 
diffusion  were  for  the  purposes  of  achieving  numerical  stability. 
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Table  4  is  a  summary  of  the  final  modeling  parameters  that  resulted  from  the 
calibration  process.  For  each  material  type,  Table  4  lists  the  following: 
Manning’s  n,  the  time-averaged  eddy  viscosity,  and  the  time-averaged  diffusion 
coefficient.  The  time-averaged  values  are  averaged  over  the  month  of  October 
1997,  but  averaging  over  any  month  in  the  simulation  period  will  yield  similar 
results  for  these  parameters.  Note  that  the  extremely  high  value  of  eddy  viscosity 
for  material  type  20  was  applied  to  ensure  stability  along  the  offshore  boundary. 

It  has  little  effect  on  the  mixing  in  the  rest  of  the  model. 


Table  4 

Calibrated  Parameters  for  Biscayne  Bay 

Material 

Type 

Description 

Manning’s  n 

Eddy 

Viscosity 

slugs/ft-sec 

Diffusion 

Coefficient 

ft2/sec 

1 

Offshore,  smooth  bed 

0.025 

40 

21 

2 

Middle  Bay,  grass  beds 

0.022 

52 

29 

3 

Card  Bank,  grass  beds 

0.045 

60 

30 

4 

South  Bay,  inlets 

0.040 

105 

38 

5 

Offshore  reefs 

0.036 

59 

18 

6 

Little  Card  Sound,  grass  beds 

0.028 

32 

23 

7 

Card  Sound,  grass  beds 

0.033 

24 

32 

8 

North  Bay,  Inlets  (Including 
Government  Cut  and  Bakers 
Haulover  Inlet) 

0.02 

84 

24 

10 

Canals 

0.035 

49 

35 

11 

Middle  Bay,  western  shoreline 

0.042 

11 

19 

12 

Barnes  Sound  and  Manatee 

Bay,  grass  beds 

0.028 

18 

12 

13 

Cutter  Bank 

0.038 

39 

30 

14  ! 

South  and  Southeast  Bay, 
southwest  of  Elliott  Key 

0.03 

39 

35 

15 

Featherbed  Bank,  Ragged  Keys, 
and  Soldier  Key 

0.089 

32 

15 

16 

North  Bay,  west  of  Virginia  Key 

0.02 

67 

31 

17 

North  Bay,  Julia  Tuttle 

Causeway  to  Bakers  Haulover 
Inlet 

0.021 

34 

21 

18 

North  Bay,  South  of  Julia  Tuttle 
Causeway 

0.021 

43 

24 

19 

Shoals,  Julia  Tuttle  Causeway  to 
Bakers  Haulover  Inlet 

0.08 

24 

37 

20 

Offshore  Boundary 

0.025 

1,148 

53 

The  Manning’s  n  values  given  here  are  typical  for  the  types  of  bottom  cover 
observed  in  Biscayne  Bay  (Chow  1959).  Of  note  are  the  extremely  high 
Manning’s  n  values  for  material  types  15  and  19.  Material  type  15  represents 
coral  and  other  subgrid-scale  roughness  features  in  the  immediate  vicinity  of 
Ragged  Keys  and  Soldier  Key,  as  well  as  a  shallow  seagrass  community  at 
Featherbed  Bank.  Material  type  19  represents  a  seagrass  community  in  very 
shallow  water.  Although  the  Manning’s  n  values  assigned  here  are  very  high, 
they  do  not  exceed  the  range  given  by  Chow. 
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Model  Verification 

After  calibration  was  complete,  a  verification  simulation  was  conducted  as 
follows:  the  model  was  run  for  the  full  10  months  of  the  simulation  period,  with 
the  final  8  months  (November  1997  through  June  1998)  used  for  hydrodynamic 
verification,  and  the  final  9  months  (October  1997  through  June  1998)  used  for 
salinity  verification.  The  model  was  run  with  a  time-step  of  0.5  hr.  An  initial 
salinity  field  was  assigned  to  the  model  based  on  field  data  collected  at  the  12 
BNP/CHL  gauges.  The  model  was  allowed  to  run  through  September  1997  to 
ensure  that  the  model  had  reached  a  state  of  dynamic  equilibrium  similar  to  that 
expected  in  the  field.  This  process  is  called  spin-up.  In  a  tidal  system,  hydraulic 
properties  propagate  at  the  speed  of  the  tidal  wave,  and  hence  reach  a  spin-up 
state  rapidly,  usually  in  a  few  tidal  cycles.  Properties  that  travel  with  the  speed 
of  a  water  particle,  such  as  salinity,  take  a  much  longer  time  to  spin  up.  Thus, 
when  salinity  is  considered,  longer  spin-up  times,  on  the  order  of  100  tidal 
cycles,  are  not  uncommon. 

Figure  22  is  a  plan  view  of  the  velocity  field  for  a  typical  ebb  tide  at 
Government  Cut.  The  contours  are  contours  of  salinity.  Note  the  offshore 
current,  which  is  induced  by  the  applied  slope  along  the  offshore  boundary  (as 
noted  in  Chapter  2).  This  current  is  responsible  for  the  observed  deflection  to  the 
north  of  the  ebb  flow  through  Government  Cut.  Figure  23  is  a  plan  view  of  the 
salinity  contours  for  April  15,  1998,  observed  at  maximum  ebb. 


Qualitative  Analysis 

Tides 

Plots  of  the  model  and  field  tides  at  the  12  BNP/CHL  CTDs  for  the  cali¬ 
bration  period  are  given  in  Plates  9  through  15,  and  the  tides  for  the  verification 
period  are  given  in  Plates  1 6  through  22.  Both  the  model  and  field  tidal  signals 
have  been  averaged  over  time,  and  the  mean  values  have  been  removed.  This 
was  done  to  more  easily  compare  the  tidal  amplitude  in  both  data  sets.  Plates  23 
through  32  contain  data  that  has  been  filtered  to  remove  the  tidal  signals.  (All 
signals  with  a  period  less  than  36  hr  were  removed).  This  was  done  to  compare 
the  effects  of  wind  on  the  water-surface  elevation  in  the  model  and  the  field 
(subtidal  effects).  The  analysis  was  conducted  for  the  months  of  October  1997, 
and  February  and  April  1998.  However,  since  the  results  were  largely  the  same 
for  both  February  and  April,  only  the  October  and  April  plots  are  given  in  this 
report.  The  analysis  was  conducted  for  CDTs  1,  3, 4,  6,  8,  9,  and  11.  In  addition, 
Plates  27  and  32  give  the  difference  between  stations  1  and  3,  and  stations  6  and 
9.  This  was  done  to  remove  the  contribution  of  the  wind  signal  imposed  at  the 
boundary,  so  the  response  of  the  model  to  wind  forcing  can  be  more  readily 
observed. 

The  plots  of  the  difference  between  stations  6  and  9  in  Plates  27  and  32 
include  the  east-west  wind  component  to  illustrate  the  dependence  of  the 
observed  response  of  the  water  surface  on  the  wind  (since  CTD6  is  essentially 
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due  west  of  CTD9,  the  east-west  component  of  the  wind  is  appropriate  for  this 
purpose). 

The  agreement  between  model  and  field  is  generally  good.  There  are  several 
instances,  however,  where  the  model  and  field  differ  in  either  tidal  amplitude  or 
the  response  of  the  signal  to  wind.  There  are  several  possible  explanations  for 
this,  including: 

a.  Errors  in  modeling  the  effects  of  wind,  including  errors  in  local  wind 
speed  specification. 

b.  Errors  in  the  field  data,  especially  drift  induced  by  meter  fouling. 

c.  Errors  in  friction  specification  in  the  model,  especially  with  respect  to  the 
potential  for  seasonal  changes  in  friction  due  to  seagrass  growth  and  die¬ 
off. 

d.  Errors  inherent  in  modeling  assumptions,  such  as  the  assumption  of 
spatially  uniform  rainfall  or  temporally  and  spatially  invariant 
atmospheric  pressure. 

The  response  of  the  water  surface  to  the  wind  is  much  greater  in  April  than  in 
October.  Therefore,  some  of  the  discrepancies  observed  in  the  April  data  were 
not  as  apparent  during  the  calibration  period  (October). 

An  overview  of  the  water-surface  elevation  data  at  each  gauge  is  given  as 
follows: 

CTD1:  The  signal  at  this  gauge  is  significantly  damped  by  friction  as  the 
tide  propagates  through  Card  Sound  and  Barnes  Sound  and  into  Manatee  Bay. 
Subtidal  signals  are  readily  apparent  in  these  data.  The  model  predicts  the  tidal 
amplitude  adequately.  Plate  28  appears  to  show  that  the  response  of  the  model  to 
the  wind  is  too  great  at  CTD1 .  However,  this  elevated  response  has  actually 
propagated  in  from  CTD3.  Plate  32  demonstrates  that  the  net  model  response  of 
the  water-surface  elevation  to  the  wind  from  CTD3  to  CTD 1  is  close  to  the  field 
response.  The  elevated  response  propagating  in  from  CTD3  is  likely  a  conse¬ 
quence  of  applying  the  tidal  data  at  CTD9  as  the  boundary  condition  across  the 
entire  offshore  boundary. 

CTD2:  The  response  at  CTD2  is  similar  to  that  observed  at  CTD1. 

CTD3:  The  tidal  signal  amplitude  in  the  model  at  this  gauge  is  consistently 
greater  than  the  amplitude  of  the  field  signal.  Additionally,  Plate  28  shows  that 
the  response  to  the  wind  is  significantly  greater  than  that  observed  in  the  field,  as 
is  mentioned  in  the  discussion  of  CTD1.  Both  of  these  discrepancies  are  likely 
consequences  of  applying  the  tidal  data  at  CTD9  as  the  boundary  condition 
across  the  entire  offshore  boundary. 

CTD4:  The  tidal  amplitude  in  the  model  and  field  are  in  good  agreement  at 
this  gauge.  Plate  28  shows  that  the  wind-induced  response  in  the  model  and  the 
field  is  in  good  agreement,  with  some  increased  response  in  the  model  during 
high-wind  events. 
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CTD5:  The  signal  in  the  field  at  this  gauge  appears  to  have  been  influenced 
by  some  type  of  local  interference.  There  is  a  tendency  for  the  field  signal  to 
truncate  the  troughs  of  the  tidal  signal.  This  may  be  the  result  of  some  error  in 
the  processing  of  the  raw  field  data.  As  of  the  publication  of  this  report,  this 
possibility  is  still  being  investigated.  Apart  from  this  aberration,  which  occurs 
throughout  much  of  the  record,  the  model  signal  at  CTD5  matches  well  with  the 
field  signal. 

CTD6:  These  data  for  this  gauge  correspond  well  with  the  field.  Of  note  is 
a  significant  deviation  between  the  model  and  field  observed  in  October. 
However,  this  deviation  is  due  to  an  instance  of  meter  fouling  which  was  left 
uncorrected  (McAdory  2002)  and  hence,  does  not  imply  modeling  error.  The 
wind-induced  response  matches  well  with  the  field.  The  field  response  observed 
toward  the  end  of  the  month  is  merely  the  meter  fouling  mentioned  previously. 
Also,  Plate  32  indicates  that  the  model  predicts  the  water-surface  slope  across  the 
bay  reasonably  well. 

CTD7:  The  agreement  between  model  and  field  is  generally  good  for  this 
gauge. 

CTD8:  The  agreement  at  this  gauge  is  similar  to  that  observed  at  CTD7. 
Plate  30  indicates  that  the  response  to  of  the  model  to  the  wind  is  also  well 
represented. 

CTD9:  The  agreement  is  quite  good  here,  which  should  be  expected  since 
the  signal  at  this  gauge  was  used  to  drive  the  model  boundary.  Plate  30  shows 
excellent  agreement  between  the  model  and  field,  indicating  that  the  decision  to 
omit  wind  stress  from  the  ocean  (to  avoid  artificially  magnifying  the  wind 
response)  was  justified. 

CTD10:  The  agreement  between  model  and  field  is  generally  good  for  this 
gauge. 

CTD11:  The  agreement  between  model  and  field  at  this  gauge  is  good. 

Plate  3 1  indicates  that  the  response  of  the  model  to  the  wind  is  also  in  good 
agreement  with  the  field. 

CTD12:  The  agreement  between  model  and  field  at  this  gauge  is  generally 
good.  The  observed  difference  in  mean  elevation  in  June  is  a  consequence  of 
meter  fouling  which  was  left  uncorrected  (McAdory  2002). 


Flows  (discharges) 

Plots  of  comparisons  between  the  model  and  the  field  flows  for  the  calibra¬ 
tion  period  are  given  in  the  Plates  33  through  45,  and  for  the  verification  period 
are  given  in  Plates  46  through  58.  In  general,  flood  tide  is  defined  as  positive 
flow,  and  ebb  tide  is  defined  as  negative  flow.  For  the  most  part,  the  agreement 
between  the  model  and  field  is  good,  especially  at  the  major  transects,  where  the 
largest  volume  exchanges  occur  (1-3,  5-8,  10-16,  20,  26).  However,  there  are 
some  significant  discrepancies.  The  flow  at  Transect  20  (Atlantic  Intracoastal 
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Waterway  (AIWW)  at  Dodge  Island  Bridge)  is  somewhat  out  of  phase  for  both 
October  and  February.  This  may  be  a  real  phase  problem,  or  it  may  be  an  error 
in  data  processing.  The  agreement  at  Transects  2  (Angelfish  Creek)  is  suspect  in 
February.  However,  the  total  volume  of  exchange  with  Biscayne  Bay  through 
Transects  2,  3  (Broad  Creek),  and  4  (Old  Rhodes  Key  to  Swan  Key)  appears  to 
be  similar  in  the  model  and  the  field. 

There  are  significant  discrepancies  observed  in  the  smaller  flows.  These 
could  be  a  consequence  of  any  of  several  factors  which  are  as  follows: 

a.  The  high  volume  of  unobserved  flow  slipping  past  along  the  boundary  of 
the  transect  (as  compared  to  the  volume  captured  by  the  transect).  This 
is  a  consequence  of  the  crudeness  of  the  shoreline  resolution  at  extremely 
small  scales. 

b.  The  low  resolution  of  the  model  for  these  flows  (often  these  transects  are 
only  one  element  wide). 

c.  The  low  velocities  of  these  flows,  resulting  in  less  dominant  advection 
and  hence,  more  erratic  flows. 


Velocities 

Plates  59  and  60  give  the  surface,  middle,  and  bottom  X-velocity  component 
for  each  of  the  ADPs,  filtered  to  remove  the  tidal  signal  (i.e.,  periods  less  than 
36  hr).  The  east- west  wind  is  included  in  the  plots  to  demonstrate  the  influence 
of  the  wind  on  the  velocities.  Although  these  plots  show  significant  differences 
in  the  response  of  each  portion  of  the  water  column  to  the  wind,  the  direction  of 
the  response  is  generally  the  same  throughout  the  water  column.  This  indicates 
that  a  depth-averaged  representation  of  the  velocities  is  a  valid  means  of  reducing 
the  data  without  losing  significant  information  about  the  behavior  of  the 
velocities. 

Note  that,  in  all  of  these  velocity  plots,  the  surface  velocity  tends  to  respond 
the  least  to  the  wind  forcing.  This  is  counterintuitive  and  could  indicate  an  error 
in  the  ADP  data,  or  in  the  reduction  of  these  data.  However,  since  no  errors  have 
been  identified  as  of  the  publication  of  this  report,  these  data  are  accepted  as  is. 

To  get  an  estimate  of  the  depth-averaged  velocity  at  each  meter  (for  purposes 
of  model  comparison),  a  weighted  average  of  the  velocity  profile  is  calculated. 
Since  the  bottom  23.6  in.  of  the  water  column  are  unmeasured,  a  logarithmic 
velocity  profile  is  assumed  in  this  region,  with  the  velocity  measured  in  the  first 
bin  (the  near  bottom  bin)  used  to  specify  the  velocity  at  the  top  of  the  profile. 
Similarly,  the  unmeasured  portion  of  the  water  column  near  the  surface  is 
assumed  to  travel  at  the  same  velocity  as  the  velocity  in  the  near  surface  bin. 

To  express  the  velocity  vector  data  set  as  a  scalar  data  set,  a  dominant  direc¬ 
tion  of  flow  was  identified  for  each  ADP.  The  X  and  Y  data  from  both  the  model 
and  field  were  used  to  generate  variance  ellipse  statistics.  These  statistics  yield 
the  magnitude  and  direction  of  the  maximum  and  minimum  variances  of  the 
velocity  data.  The  dominant  direction  of  flow,  maximum  standard  deviation 
(taken  from  the  maximum  variance),  minimum  standard  deviation,  and  the  ratio 
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of  the  standard  deviations  (maximum  to  minimum)  for  each  depth-averaged  ADP 
data  record  are  given  in  Table  5.  The  ratio  of  the  standard  deviations  indicates 
the  degree  of  uniformity  in  the  direction  of  the  velocities;  a  high  ratio  indicates 
that  the  velocities  oscillate  along  a  nearly  constant  path,  whereas  a  low  ratio 
indicates  that  the  velocities  oscillate  along  multiple  paths.  Hence,  high  ratios 
indicate  confined  flow,  and  low  ratios  indicate  unconfined  flow. 


Table  5 

Variance  Ellipse  Statistics 

ADP 

Dominant  Direction  of  Flow 
(degrees  ccw  from  the 
positive  x-axis) 

Maximum 

Standard 

Deviation 

ft/sec 

Minimum 

Standard 

Deviation 

ft/sec 

Standard 

Deviation 

Ratio 

ADP1  MODEL 

42.32 

0.29 

0.172 

1.69  ‘ 

ADP1  FIELD 

2.21 

0.352 

0.198 

1.78 

ADP2  MODEL 

82.8 

0.27 

0.061 

4.44 

ADP2  FIELD 

68.55 

0.319 

0.235 

1.36 

ADP3  MODEL 

23.48 

0.696 

0.145 

4.81 

ADP3  FIELD 

18.92 

0.5 

0.206 

2.43 

ADP5  MODEL 

48.03 

0.801 

0.023 

35.02 

ADP5  FIELD 

39.71 

0.676 

0.163 

4.14 

The  meters  were  active  throughout  the  entire  period  of  record.  ADP  1,  how¬ 
ever,  was  inactive,  from  September  through  December  1 997.  Therefore,  there 
are  no  data  at  ADP1  for  the  calibration  period. 

The  velocities  at  ADP’s  2,  3,  and  5  are  plotted  over  the  calibration  period  in 
Plates  61  and  62.  The  plots  of  the  velocities  over  the  verification  period  are 
plotted  in  Plates  63  and  64.  For  each  plot,  the  dominant  directions  taken  from 
Table  5  are  used  to  determine  a  flow  axis.  Each  velocity  measurement  is  then 
projected  onto  the  flow  axis,  and  the  sign  of  the  projection  is  used  to  determine 
whether  the  flow  is  positive  or  negative.  Plates  65  through  71  contain  data  which 
have  been  filtered  to  remove  the  tidal  signals  (all  signals  with  a  period  less  than 
36  hr  were  removed).  This  was  done  to  compare  the  effects  of  wind  on  the 
velocities  in  the  model  and  the  field.  The  analysis  was  conducted  for  the  months 
of  October  1997  and  February  and  April  1998.  However,  since  the  results  were 
largely  the  same  for  both  February  and  April,  only  the  October  and  April  plots 
are  given  in  this  report.  For  these  plots,  the  velocity  components  were  plotted 
(i.e.,  X  and  Y  velocity). 

Some  comments  concerning  the  velocities  at  each  ADP  meter  are  given 
below. 

ADP1:  The  model  appears  to  simulate  the  tidal  velocity  amplitudes  with 
reasonable  accuracy.  Plate  68  indicates  that  the  model  is  in  good  agreement  with 
the  field  response  to  the  wind.  Note  that  Table  5  gives  a  dominant  direction  of 
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flow  for  the  field  that  has  a  greater  east-west  component  than  the  model.  Since 
the  wind-induced  currents  are  well  represented  in  the  model,  this  difference  in 
the  dominant  direction  of  flow  is  likely  a  difference  in  the  direction  of  the  tidal 
excursion.  This  difference  may  be  due  to  a  local  subgrid-scale  geometric  feature 
in  the  field  that  is  not  resolved  in  the  model. 

ADP2:  The  model  appears  to  simulate  the  tidal  velocity  amplitudes  with 
reasonable  accuracy.  However,  Plates  65  and  69  show  that  the  wind  response  in 
the  model  and  the  field  are  in  the  opposite  directions,  and  of  drastically  different 
magnitudes.  A  possible  explanation  for  this  can  be  inferred  from  Figure  24.  This 
is  a  plot  of  the  residual  velocities  in  the  Central  Bay,  averaged  over  April  15  and 
16, 1998.  Note  that,  in  the  vicinity  of  ADP2,  the  model  predicts  a  residual 
recirculation  current.  If  the  location  of  the  meter  in  the  field  is  different  than  the 
model,  or  if  the  local  geometry  and/or  mixing  as  represented  in  the  model  are 
slightly  different,  then  the  residual  current  observed  in  the  model  could  reverse 
direction  to  correspond  with  the  field. 

ADP3:  The  amplitude  of  the  tidal  velocity  is  greater  in  the  model  than  in  the 
field.  This  is  surprising  because  the  measured  flow  through  the  Safety  Valve, 
which  is  located  just  east  of  ADP3,  correlates  extremely  well  with  that  given  in 
the  field  (see  Transect  8,  Plates  36  and  49).  Therefore,  this  effect  is  likely  local. 
The  wind-induced  current  at  ADP3  corresponds  well  with  that  observed  in  the 
field  (Plates  66  and  70). 

ADP5:  The  tidal  velocities  in  the  model  at  this  meter  correspond  extremely 
well  to  those  observed  in  the  field.  The  wind-induced  currents  are  also  in  good 
agreement  with  the  wind-induced  currents  in  the  field  (Plates  67  and  71). 


Salinities 

The  salinity  plots  at  the  12  BNP/CHL  CTDs  are  given  for  the  verification 
period  in  Plates  72  through  77.  These  plates  also  contain  selected  plots  of 
DERM  Bay  Run  salinity  data  for  Bay  Run  gauges  with  are  in  close  proximity  to 
the  BNP/CHL  CTDs. 

In  general,  the  model  and  field  are  in  good  agreement  with  respect  to  salinity. 
However,  where  discrepancies  exist,  they  are  likely  the  result  of  any  of  several 
factors.  These  factors  are  (in  order  of  likely  significance): 

a.  Inaccurate  values  of  local  wind  strength  applied  in  the  model. 

b.  Inaccurate  or  erroneous  estimates  of  the  inflow  rates  at  the  coastal 
structures. 

c.  Insufficient  or  excessive  currents  in  the  offshore  to  sweep  away  low- 
salinity  water. 

d.  Errors  in  the  estimation  of  the  offshore  salinity  boundary  condition. 

e.  Locally  stratified  conditions  in  the  bay  (since  the  gauges  are  near  the 
bed,  locally  stratified  conditions,  if  and  when  they  exist,  would  tend  to 
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make  the  observed  values  higher  than  the  depth  averaged  value  given  in 
the  model). 

f  Errors  in  the  estimation  of  the  evaporation  rate. 

g.  The  lack  of  any  tidal  wetlands  resolved  in  the  grid,  which  would  supply 
an  additional  evaporation  sink  for  fresh  water. 

A  discussion  of  each  gauge  follows: 

CTD1:  The  model  reasonably  predicts  the  field  measurement  for  much  of 
the  verification  period.  However,  the  model  exceeds  the  field  value  by  as  much 
as  3  to  4  from  December  through  March.  This  may  be  due  to  the  excessive  wind 
response  in  the  tide  propagating  in  from  the  offshore  boundary,  as  observed  at 
CTD3  (see  “Tides”  section).  Note  that  a  large  storm  event  recorded  in  early 
February  is  manifest  in  both  the  model  and  in  the  field.  The  model  tends  to 
underpredict  the  amplitude  of  the  salinity  signal,  possibly  indicating  excessive 
turbulent  mixing  in  the  model.  However,  Table  4  indicates  an  average  diffusion 
coefficient  for  the  material  type  in  Manatee  Bay  (material  type  12)  of  12  fr/s, 
which  is  relatively  small. 

CTD2:  This  gauge  is  influenced  by  flow  at  S-20.  There  is  a  large  flow 
event  recorded  at  S-20  in  September.  The  flow  rate  recorded  for  this  event 
reaches  900  cfs,  which  exceeds  the  Special  Project  Flood  for  this  structure 
(750  cfs)  (Swain  et  al.  1997).  This  may  indicate  an  error  in  the  estimate  of  this 
flow.  The  model  tends  to  overestimate  the  salinity  in  February  and  March, 
possibly  because  of  the  influence  of  excessive  wind  response  in  the  tide  (see  the 
discussion  of  CTD1).  In  May,  the  salinity  in  the  field  climbs  to  two  (2)  to  three 
(3)  parts  per  thousand  (ppt)  greater  than  salinity  than  the  model  predicts.  This 
increase  in  salinity  is  also  observed  at  CTD4,  CTD5,  and  CTD6.  This  rapid 
increase  in  salinity  in  the  bay  is  the  result  of  minimal  inflow  from  the  coastal 
structures,  combined  with  a  prevailing  onshore  wind.  The  model  predicts  the 
trend  well,  but  tends  to  slightly  underpredict  the  rate  of  increase  in  the  salinity. 

CTD3:  The  model  is  within  1  to  2  ppt  of  the  field  for  much  of  the  record, 
with  a  3-  to  4-ppt  difference  through  October  and  November.  Note  that  there  is 
some  slight  stratification  observed  in  the  Seabird  data  during  this  period,  which 
may  explain  some  of  the  discrepancy. 

CTD4:  These  data  at  this  gauge  are  plotted  together  with  data  from  two  of 
the  DERM  Bay  Run  stations:  BB41,  located  about  1.5  miles  south  of  CTD4,  and 
BB38,  located  about  2.5  miles  northeast  of  CTD4.  The  model  and  field  are  in 
good  agreement  for  most  of  the  record,  with  a  departure  of  2  to  3  ppt  occurring  in 
April  and  May  (see  discussion  of  CTD2).  The  DERM  Bay  Run  data  show  some 
intermittent  stratification  in  the  region.  Of  note  is  stratification  observed  at  BB41 
in  early  February,  which  corresponds  to  a  major  flood  event.  This  event  is  also 
observed  in  the  model  data. 

CTD5:  Although  the  magnitude  of  the  salinity  for  the  model  and  field  are 
often  different,  the  pattern  of  the  salinity  response  is  well  represented  in  the 
model.  This  response  is  largely  due  to  the  influence  of  wind  and  the  inflows  at 
the  nearby  canals.  It  is  noteworthy  that  this  region  of  the  bay  exhibits  high 
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horizontal  salinity  gradients.  Therefore,  the  value  of  salinity  recorded  at  nearby 
locations  can  be  significantly  different.  Also,  the  observed  error  in  the  predicted 
inflow  for  the  nearby  coastal  structures  (i.e.,  S-20F,  S-20G,  S-21,  and  S-21A)  is 
relatively  high  (Table  3).  Therefore,  errors  in  the  salinity  can  be  partially 
attributed  to  errors  in  the  estimate  of  the  inflows  in  the  canals. 

These  data  at  this  gauge  are  plotted  together  with  data  from  the  DERM  Bay 
Run  station  MI01,  which  is  located  in  the  mouth  of  Military  Canal.  It  is  note¬ 
worthy  that  S-20G  is  rarely  flowing  though  the  simulation  period  (Plate  2). 
Therefore,  the  high  stratification  observed  at  MI01  is  likely  the  result  of  flows 
through  the  other  nearby  canals.  Although  the  Seabird  data  indicate  little  strati¬ 
fication,  the  shallow  depth  at  this  gauge  (average  depth  is  2.9  ft)  makes  the 
measurement  of  a  salinity  profile  with  the  seabird  device  difficult  (the  Seabird 
measurement  excludes  the  upper  3  ft  of  the  water  column).  Hence,  the  Seabird 
measurement  at  this  gauge  is  not  helpful  in  determining  the  degree  of  local 
stratification  at  CTD5. 

CTD6:  These  data  at  this  gauge  are  plotted  together  with  data  from  the 
DERM  Bay  Run  station  BB39A,  which  is  located  within  a  few  hundred  feet  of 
CTD6.  The  salinity  at  CTD6  is  generally  in  good  agreement  with  the  field,  both 
with  respect  to  average  salinity  and  salinity  amplitude.  However,  the  model 
tends  to  overpredict  the  salinity  amplitude  during  January  through  March.  There 
is  some  local  stratification  observed  at  BB39A.  This  gauge  is  also  influenced  by 
inflows  from  S-20F,  S-20G,  S21-A,  and  S-21,  although  to  a  lesser  extent  than  is 
CTD5. 

CTD7:  This  gauge  is  located  within  2  miles  of  S-123.  It  is  plotted  together 
with  data  from  the  DERM  Bay  Run  stations  CD01A,  which  is  located  within  a 
few  hundred  feet  of  CTD7,  and  BB36,  which  is  about  3  miles  east  of  CTD7.  In 
general,  the  agreement  between  the  model  and  field  is  good,  except  that  the 
model  tends  to  overpredict  salinity  through  December  and  January  by  about  3  to 
4  ppt.  The  effects  of  the  February  storm  event  are  manifest  in  both  the  model 
and  the  field.  Note  that,  toward  the  end  of  the  period  of  record,  CD01A  records 
persistent  high  stratification,  whereas  the  Seabird  gauge  records  none.  Since 
CD01A  is  located  nearby  CTD7,  this  points  to  the  possibility  that  the  Seabird 
gauge  is  not  recording  local  stratification.  If  this  is  true,  it  may  be  because  the 
freshwater  lens  is  too  close  to  the  surface  to  be  measured  by  the  Seabird  gauge. 
Note  also  the  large  salinity  difference  between  CD01A  and  BB36,  which  are 
separated  by  about  3  miles.  These  differences  demonstrate  the  high  horizontal 
salinity  gradients  near  the  western  shoreline.  This  is  especially  evident  from 
April  through  June. 

CTD8:  This  gauge  is  located  in  the  mouth  of  Canal  2,  which  is  controlled  by 
S-22  (Plate  4).  This  gauge  is  plotted  together  with  data  from  Bay  Run  station 
BB34,  which  is  located  within  a  few  hundred  feet  of  CTD8.  The  model  and  field 
are  in  close  agreement  for  the  entire  period  of  record.  It  is  noteworthy  that, 
although  the  Seabird  data  evidences  little  or  no  stratification,  the  BB34  data 
record  persistent  stratification  throughout  the  period  of  record.  This  is  further 
evidence  to  support  the  hypothesis  that  the  Seabird  gauge  is  not  measuring  the 
fresher  water  at  the  surface.  There  is  particularly  strong  stratification  observed  at 
the  time  of  the  rainfall  event  in  February. 


Chapter  3  Model  Calibration  and  Verification 


29 


CTD9:  This  gauge  agrees  reasonably  well  with  the  field,  although  the 
amplitude  of  the  signal  is  generally  smaller  than  that  observed  in  the  field.  The 
salinity  at  this  gauge  is  largely  a  reflection  of  the  value  specified  at  the  offshore 
boundary.  Therefore,  the  close  agreement  at  this  gauge  tends  to  support  the 
validity  of  the  choice  of  the  salinity  boundary  condition  for  the  offshore 
boundary. 

CTD10:  This  gauge  is  plotted  together  with  several  nearby  DERM  Bay  Run 
stations.  CTD10  is  located  near  the  mouth  of  the  Miami  River.  The  model 
agrees  well  with  the  field,  in  both  mean  salinity  and  salinity  tidal  amplitude. 
DERM  Bay  Run  station  MR01  is  located  in  the  mouth  of  the  Miami  River,  and 
hence,  records  significant  stratification  throughout  the  period  of  record.  There  is 
evidence  of  meter  fouling  in  the  field  data  toward  the  end  of  June. 

CTD11:  This  gauge  agrees  well  with  the  field.  There  is  evidence  of  meter 
fouling  in  the  field  data  toward  the  end  of  June. 

CTD12:  Although  the  agreement  between  model  and  field  at  CTD12  is 
generally  good,  there  is  a  significant  drop  in  model  salinity  observed  in  early 
December  and  in  early  February.  There  is  a  high  flow  record  at  S-28  that  occurs 
in  these  time  frames,  and  hence,  errors  in  the  estimate  of  flow  at  this  location, 
and  this  could  influence  CTD12.  Also,  there  is  relatively  high  stratification 
observed  at  the  local  DERM  Bay  Run  stations  and  in  the  Seabird  data  at  these 
times,  which  could  contribute  to  the  discrepancy. 


Quantitative  Analysis 

Statistical  analysis  of  the  model  and  field  data  is  given  with  four  statistics: 
mean  error  (ME),  mean  absolute  error  (MAE),  root  mean  square  error  (RMS), 
and  index  of  agreement  (d).  The  index  of  agreement  is  a  descriptive  measure  of 
predictive  success.  It  gives  a  value  between  0  and  1,  with  1  indicating  perfect 
agreement.  Since  it  is  a  relative  measure  of  modeling  success,  and  it  is  bounded 
(between  0  and  1 ),  it  is  also  useful  for  making  cross-comparisons  between 
models  (Willmott  1982). 

Letting  O  and  P  be  the  observed  and  predicted  values  of  a  quantity  N,  the 
following  expressions  result: 

ME=^f(0„-P„)  (1) 

(2) 

/V  n~] 


RMS  = 


liP'-P. 


) 


(3) 
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The  agreement  between  the  model  and  field  is  very  high.  This  is  evident  in 
the  values  given  for  the  index  of  agreement,  which  are  very  close  to  1 .  The  best 
agreement  is  at  CTD9  and  Virginia  Key,  which  is  not  surprising  since  data  at 
CTD9  were  used  to  drive  the  offshore  boundary,  and  the  NOS  gauge  is  located 
nearby.  The  worst  agreement  occurs  at  CTD  Gauges  1  through  3  and  Gauge  5. 
The  errors  at  Gauges  1  through  3  are  likely  propagated  in  from  the  offshore 
boundary,  which  was  specified  according  to  the  tide  recorded  at  CTD9,  located 
further  to  the  north.  Note  that  the  error  shown  for  CTD5  may  be  exaggerated 
because  of  the  truncation  of  the  troughs  in  the  field  record  that  was  discussed  in 
the  “Qualitative  Analysis”  section. 

A  harmonic  analysis  was  performed  on  the  water-surface  elevation  data.  The 
analysis  was  done  over  the  month  of  May  1998  because  there  exists  field  data  at 
all  the  gauges  for  this  month.  The  tidal  signal  was  filtered  to  remove  all  signals 
with  a  period  less  than  3  hr  and  greater  than  35  hr.  The  signal  was  then  analyzed 
for  six  constituents:  M2,  S2,  N2,  K2,  Kl,  and  01.  The  results  are  given  in 
Plates  78  through  85.  Plates  78  through  84  compare  the  tidal  constituent  ampli¬ 
tudes.  These  confirm  that  the  M2  constituent  is  the  dominant  component  of  the 
tide.  Plate  85  gives  the  phase  lag  of  the  constituents,  relative  to  CTD9.  For 
convenience,  the  average  phase  lag  is  calculated  for  the  semidiurnal  components 
(M2,  S2,  N2,  K2)  and  the  diurnal  components  (Kl,  01).  The  agreement  in 
amplitude  and  phase  is  generally  good.  The  worst  agreement  appears  at  CTD3. 
This  may  result  from  the  fact  that  the  tide  data  at  CTD9  were  used  as  a  model 
boundary  condition  across  the  entire  offshore  boundary.  Note  that,  although  the 
phase  lag  is  especially  poor  for  the  diurnal  components  at  CTD3,  there  is  very 
little  tidal  energy  associated  with  the  diurnal  components. 


Velocities 

A  statistical  analysis  of  the  velocity  data  for  the  calibration  period  is  given  in 
Table  8,  and  Table  9  contains  the  same  analysis  for  the  verification  period. 

These  tables  contain  the  same  statistics  that  are  given  in  Tables  6  and  7  for  the 
water-surface  elevations.  The  statistics  bear  out  the  qualitative  analysis:  the  best 
agreement  between  model  and  field  is  at  ADP5,  the  worst  is  at  ADP2. 


Salinities 

A  statistical  summary  of  the  salinity  data  at  the  BNP/CHL  gauges  is  given  in 
Table  10.  The  summary  contains  the  same  statistics  as  those  given  for  the  water- 
surface  elevations  and  velocities.  Although  the  mean  difference  in  salinity  at 


Table  8 

Statistics  for  Velocity  Comparisons:  Calibration 


ADP 

ME,  ft/sec 

MAE,  ft/sec 

RMS,  ft/sec 

d 

1 

No  Data 

No  Data 

No  Data 

No  Data 

2 

-0.026 

0.185 

0.253 

0.842 

3 

-0.016 

0.269 

0.334 

0.924 

5 

-0.028 

0.138 

0.184 

0.984 

32 


Chapter  3  Model  Calibration  and  Verification 


Table  9 

Statistics  for  Velocity  Com 

parisons:  Verification 

ADP 

ME  (ft/sec) 

MAE  (ft/sec) 

RMS  (ft/sec) 

d 

1 

-0.006 

0.179 

0.258 

0.879 

2 

-0.041 

0.275 

0.35 

0.714 

3 

-0.032 

0.208 

0.272 

0.954 

5 

-0.036 

0.17 

0.246 

0.973 

Table  10 
Statistics  for 

Salinity  Comparison 

CTD  Gauge 

ME,  ppt 

MAE,  ppt 

RMS,  ppt 

d 

1 

-1.154 

1.528 

1.962 

0.657  1 

2 

0.119 

1.927 

2.304 

0.594 

3 

1.607 

2.225 

2.922 

0.496 

4 

0.177 

1.517 

1.858 

0.812 

5 

3.643 

4.166 

4.967 

0.886 

6 

-0.396  ] 

2.271 

2.922 

0.931 

7 

-1.105 

1.963 

2.562 

0.867 

8 

0.612 

1.264 

1.677 

0.94 

9 

-0.154 

0.847 

1.398 

0.718 

10 

-1.2 

2.345 

3.784 

0.639 

11 

0.406 

1.103 

2.232 

0.794 

12 

1.806 

2.496 

3.624 

0.501 

CTD5  is  relatively  high,  this  gauge  also  exhibits  a  high  index  of  agreement.  This 
is  because  of  the  ability  of  the  model  to  predict  the  general  behavior  of  the 
salinity  in  the  field.  CTD  Gauges  1  through  3  yield  a  low  index  of  agreement. 
This  is  a  result  of  the  low  variability  of  the  salinity  signal  over  time,  combined 
with  the  tendency  for  the  model  to  overpredict  the  salinity  in  the  middle  of  the 
simulation.  The  gauges  in  the  North  Bay  show  higher  values  of  mean  absolute 
error  than  of  mean  error,  indicating  the  dominance  of  salinity  amplitude  differ¬ 
ences,  or  the  lack  of  a  persistent  trend  in  overprediction  or  underprediction. 


Chapter  3  Mode!  Calibration  and  Verification 


33 


4  Salinity  Sensitivity  Tests 


Five  tests  were  conducted  to  determine  the  sensitivity  of  the  model  salinity  to 
variations  in  the  boundary  conditions.  The  tests  were  designed  to  reflect  a  con¬ 
servative  estimate  of  the  degree  of  uncertainty  in  the  model  boundary  forcings. 
The  tests  are  as  follows: 

a.  Test  A:  Sensitivity  to  the  imposed  offshore  current  was  investigated  by 
eliminating  the  offshore  boundary  slope. 

b.  Test  B:  Sensitivity  to  the  evaporation  rate  was  investigated  by  increasing 
the  evaporation  rate  by  a  factor  of  1 .5. 

c.  Test  C:  Sensitivity  to  errors  in  the  estimates  of  the  flows  at  the  structures 
was  investigated  by  reducing  the  amount  of  flow  at  each  structure  by 
either 

(1)  One-half  the  error  given  in  Table  3,  or 

(2)  10  percent,  whichever  is  larger. 

d.  Test  D :  Sensitivity  to  the  offshore  salinity  boundary  condition  was 
investigated  by  increasing  the  specified  boundary  condition  by  1  ppt. 

e.  Test  E:  Sensitivity  to  the  wind  speed  (and  hence  wind  stress)  was 
investigated  by  increasing  the  wind  speed  by  a  factor  of  1 .2. 

The  results  of  these  sensitivity  tests  are  given  in  Figure  25.  The  figure  shows 
the  mean  absolute  value  of  the  change  in  salinity  at  each  CTD  gauge  because  of 
each  sensitivity  test.  The  most  significant  impact  results  from  the  increase  in 
wind  speed  (Test  E).  This  impact  is  most  evident  along  the  western  shoreline  of 
Biscayne  Bay,  at  CTDs  5,  6,  7,  and  8.  The  reduction  in  the  inflow  rates  at  the 
coastal  structures  (Test  C)  has  a  significant  impact  on  the  salinity  at  CTDs  5  and 
6,  where  the  flows  at  structures  S-20F,  S20-G,  S21,  and  S-21A  have  a  dramatic 
influence.  The  increased  evaporation  (Test  B)  has  a  significant  impact  in  South 
and  Central  Biscayne  Bay.  The  elimination  of  the  offshore  current  (Test  A)  has 
the  greatest  impact  at  CTD3.  The  increased  offshore  salinity  (Test  D)  has  the 
greatest  impact  in  the  North  Bay. 

Table  1 1  gives  the  vector  sum  of  the  mean  absolute  differences  in  salinity  for 
the  five  sensitivity  tests.  This  is  given  together  with  the  mean  absolute  error  in 
salinity  at  each  of  the  gauges  for  the  verification  run.  If  the  difference  between 
these  two  values  is  positive  (i.e.,  if  the  vector  sum  of  sensitivity  differences  is 
greater  than  the  mean  absolute  error),  we  can  assume  that  the  error  in  the  model 
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Table  11 

Comparison  of  Sensitivity  to  Baseline  Run 

CTD  Gauge 

Vector  Sum  of 
Sensitivity 

Differences 

MAE  for  Baseline 
Model  Run 

Difference 

1 

1.43 

1.53 

-0.10 

2 

1.55 

1.93 

-0.37 

3 

1.53 

2.23 

-0.70 

4 

1.71 

1.52 

0.19 

5 

3.48 

4.17 

-0.68 

6 

2.94 

2.27 

0.67 

7 

1.91 

1.96 

-0.06 

8 

1.79 

1.26 

0.53 

9 

1.14 

0.85 

0.29 

10 

1.45 

2.35 

-0.90 

11 

1.45 

1.10 

0.34 

12 

1.27 

2.50 

-1.23 

is  within  the  bounds  of  the  error  in  the  applied  boundary  conditions.  Table  1 1 
shows  several  negative  differences,  which  are  given  in  bold  type.  However,  all  of 
these  differences  (except  for  CTD12)  are  less  than  lppt,  and  given  the  fact  that 
the  sensitivity  tests  were  conducted  with  conservative  estimates  of  the  error  in  the 
applied  boundary  conditions,  these  negative  differences  do  not  indicate  conclu¬ 
sively  that  the  model  is  performing  inadequately  at  any  of  the  gauges.  The  per¬ 
formance  of  the  model  at  CTD12  is  of  concern,  but  it  may  be  that  this  results 
from  errors  in  the  estimate  of  the  inflow  at  S-28  (p  25),  and  that  applying  the 
conservative  estimate  of  one-half  the  measured  error  in  the  inflow  may  have  been 
inadequate. 
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Hydrodynamics 

Tides:  The  model  successfully  reproduces  the  observed  tidal  signal  in  the 
Central  and  North  Bay,  with  respect  to  both  the  tidal  and  subtidal  (i.e.,  wind- 
driven)  response.  The  tidal  signal  is  well  represented  in  the  South  Bay,  but  the 
wind-driven  response  is  excessive  in  the  model.  This  appears  to  result  from  the 
application  of  the  measured  tidal  signal  at  CTD9  as  the  offshore  boundary  condi¬ 
tion.  Although  this  boundary  specification  is  valid  in  the  North  and  Central  Bay, 
it  propagates  an  exaggerated  wind  response  into  the  South  Bay. 

Flows:  The  model  adequately  reproduces  the  observed  flows  across  the 
ADCP  transects.  Although  there  are  some  discrepancies,  most  of  these  are  across 
transects  with  relatively  small  flow  rates. 

Velocities:  The  model  successfully  reproduces  the  observed  velocity  signal, 
both  tidal  and  subtidal,  at  ADPs  3  and  5.  The  subtidal  signal  is  well  represented 
at  ADP1,  but  the  dominant  direction  of  the  tidal  excursion  in  the  model  is  more 
northerly  than  the  field,  by  about  40  deg.  The  tidal  signal  is  adequately  repre¬ 
sented  at  ADP2,  but  the  subtidal  velocities  are  in  the  opposite  directions  in  the 
model  and  the  field.  This  is  likely  a  consequence  of  the  fact  that  ADP2  is  located 
in  a  residual  return  current  in  the  model,  and  the  local  configuration  of  this 
current  could  be  different  in  the  field. 

Also,  subtidal  analysis  of  the  surface,  middepth,  and  bottom  velocities  for 
ADPs  1,  2,  3,  and  5  reveals  that,  although  the  response  to  the  wind  is  different  at 
different  depths,  the  entire  water  column  generally  responds  in  the  same  direc¬ 
tion.  This  serves  as  further  justification  for  the  validity  of  depth  averaging,  and 
hence,  the  use  of  a  2-D  depth-averaged  model  for  Biscayne  Bay. 


Salinity 

The  model  adequately  reproduces  the  observed  field  salinity.  There  are  some 
significant  discrepancies  between  the  model  and  the  field,  but  sensitivity  tests 
conducted  to  determine  the  sensitivity  of  the  model  to  applied  boundary  condi¬ 
tions  tend  to  indicate  that  the  error  in  the  model  is  (a)  within  the  bounds  of  the 
uncertainty  in  the  boundary  conditions,  or  (b)  exceeds  these  bounds,  but  by  an 
insignificant  amount,  especially  considering  that  the  bounds  of  the  boundary 
uncertainty  were  chosen  conservatively. 
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Overall  Model  Assessment 


The  model  is  sufficiently  well-calibrated  for  use  in  conducting  base  versus 
plan  simulations  to  assess  the  impact  of  various  freshwater  distribution  scenarios 
on  the  salinities  in  Biscayne  Bay. 


Model  Improvements 

Some  suggested  improvements  include: 

a.  A  thorough  investigation  of  the  inflow  rates  at  the  canals. 

b.  The  implementation  of  an  ocean  boundary  condition  generated  by  a 
coastal  ocean  model  (such  as  ADCIRC)  with  the  effects  of  wind 
included. 

c.  A  field  investigation  of  the  average  wind  speed  and  direction  in  the  vari¬ 
ous  parts  of  Biscayne  Bay.  This  might  consist  of  a  program  of  synoptic 
wind  data  collection  in  several  different  locations,  using  the  same  type  of 
equipment  mounted  at  the  same  elevation  above  the  water  surface. 

d.  The  resolution  of  coastal  wetlands  into  the  grid. 


Model  Computational  Performance 

The  model  was  run  on  two  separate  systems:  a  DEC-Alpha  workstation  with 
four  processors  running  in  parallel,  and  the  Silicon  Graphics  03K  system,  using 
eight  processors  running  in  parallel.  For  a  10-month  simulation  with  one-half¬ 
hour  time-steps,  the  Dec- Alpha  required  approximately  6.5  days  to  complete  the 
run,  and  the  03K  required  approximately  1.2  days  to  complete  the  run. 
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Figure  1 .  Biscayne  Bay  location  maps 
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Figure  2.  Feature  map  of  Biscayne  Bay 


Figure  3.  The  FE  mesh 


Figure  5.  The  middle  mesh,  showing  Middle  Bay 


Figure  7.  Material  types  for  Biscayne  Bay  (as  described  in  Table  1) 


Figure  8.  Location  of  salinity  control  structures  and  offshore  salinity  gauges  for 
Biscayne  Bay 
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Figure  10.  Wind  speed  at  Convoy  Point 
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Figure  1 1 .  Wind  multiplication  factors  for  Biscayne  Bay 
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Figure  13.  Evaporation  at  Convoy  Point 
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Figure  14.  Flow  and  salinity  at  Jewfish  Creek 
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Figure  15.  Application  locations  for  groundwater  inflow 
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Figure  16.  Groundwater  inflow  to  Biscayne  Bay 


Figure  17.  Instrument  location  map 


Figure  18.  Instrument  location  map  (TR  =  Transect) 


Figure  19.  DERM/BNP  stations  in  North  Bay  (Courtesy  Tim  McIntosh,  DERM.) 


Figure  20.  DERM/BNP  stations  in  Middle  Bay  (Courtesy  Tim  McIntosh,  DERM) 


Figure  21.  DERM/BNP  stations  in  South  Bay  (Courtesy  Tim  McIntosh,  DERM) 


Figure  23.  Salinity  contours  for  April  15,  1998,  at  maximum  ebb 
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Figure  25.  Absolute  difference  in  salinity  resulting  from  sensitivity  tests 
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Appendix  A 
TABS-MDS  Theoretical 
Development 


This  appendix  contains  a  partial  development  of  the  equations  used  in  TABS- 
MDS  (Multi-Dimensional  Sediment),  as  well  as  a  description  of  some  of  the 
empirical  and  theoretical  expressions  used  to  simulate  various  flow  phenomena. 
The  development  of  the  equations  includes  the  development  of  equations  for 
three-dimensional  (3-D),  two-dimensional  (2-D)  laterally  averaged,  2-D  verti¬ 
cally  averaged,  and  one-dimensional  (1-D)  elements.  Although  only  2-D  verti¬ 
cally  averaged  elements  are  used  in  the  present  study,  the  form  of  the  equations 
results  from  their  reduction  from  the  3-D  form,  and  hence  an  adequate  explana¬ 
tion  of  the  2-D  equations  requires  the  inclusion  of  some  discussion  of  the  3-D 
equations. 


TABS-MDS  Introduction 

TABS-MDS  is  a  finite  element,  hydrodynamic  model.  It  is  based  on 
RMA10,  a  model  written  by  Ian  King,  Resource  Management  Associates.1  It  is 
capable  of  modeling  turbulent,  subcritical  flows  using  1-D,  2-D,  and/or  3-D 
elements.  It  is  also  capable  of  modeling  constituent  transport.  This  includes 
modeling  salinity,  temperature,  and/or  fine-grained  sediment.  The  model  is 
capable  of  coupling  the  spatial  density  variation  induced  by  concentration  gradi¬ 
ents  in  the  constituent  field  to  the  hydrodynamic  calculations.  This  enables  the 
model  to  simulate  phenomena  such  as  saline  wedges  in  estuaries.  The  model  has 
features  that  permit  the  simulation  of  intermittently  wetted  regions  of  the  domain, 
such  as  coastal  wetlands. 


1  All  references  cited  in  this  appendix  are  listed  in  the  References  section  at  the  end  of  the  main 
text. 
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A1 


TABS-MDS  Theoretical  Development 

3-D  equations 

There  are  six  unknowns  ( u ,  v,  w,  h,  s,  p).  Therefore,  six  equations  are 
required. 

The  Navier-Stokes  equations  (i.e.,  conservation  of  fluid  momentum) 
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The  Volume  Continuity  Equation 

du  +  dy  dvr  __ 
dx  dy  dz 
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The  Advection-Diffusion  Equation 
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(5) 


The  Equation  of  State 

p  =  F(s,/)  ^ 

where 

x  =  applied  forces  (e.g.,  wind  stress,  bed  shear  stress,  Coriolis  force) 

0j  =  salinity  source/sink  term 

Now  reduce  the  number  of  unknowns  requiring  a  simultaneous  solution  from 
6  to  3. 

Assuming  that  the  influence  of  vertical  momentum  on  the  system  is  small 
and  may  be  neglected,  Equation  3  reduces  to  the  following  equation: 


az 


(7) 


Equation  7  is  a  statement  that  the  vertical  pressure  distribution  is  hydrostatic. 

Equation  4  may  then  be  integrated  in  the  vertical  direction  to  yield  the 
following  equation: 


\  °+rhdw, 

ch i  =  -  j  -dr\  =  -ws+wb 


where 

ws  =  vertical  velocity  at  the  water  surface 
Wb  -  vertical  velocity  at  the  bed 

The  surface  velocity  can  be  expressed  as  follows: 


(8) 
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(9) 


d(zb  +  h) 
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ox 
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dy  dt 


Similarly,  the  bed  velocity  can  be  expressed  as: 


wb 


=  u. 


dx 


■  +  vk 


dzb  ,  dzh 
dy  dt 


(10) 


where 

us,  v*  =  surface  horizontal  velocity  components 
Mi,  Vi  =  near  bed  horizontal  velocity  components 
Zb  =  bed  elevation 

Note  that  by  replacing  Equations  3  and  4  with  6  and  8,  we  recast  the  equa¬ 
tions  such  that  w  is  present  only  in  the  horizontal  momentum  equations  and  the 
advection  diffusion  equation.  It  can  now  be  solved  in  a  separate  decoupled 
calculation  using  the  original  form  of  the  continuity  equation  (Equation  4).  This 
is  done  by  taking  the  derivative  of  equation  4  with  respect  to  z  and  solving  for  w, 
applying  ws  and  wb  as  boundary  conditions. 

You  can  further  eliminate  p  from  the  list  of  unknowns  requiring  a  simul¬ 
taneous  solution  by  solving  the  equation  of  state  (Equation  6)  in  a  decoupled 
step. 

Thus,  you  are  left  with  four  equations  (1,  2,  8,  and  5)  and  four  unknowns  (m, 
v,  h  and  s)  to  be  solved  simultaneously.  In  practice,  however,  the  solution  is 
broken  up  into  two  steps.  First  the  velocities  and  depth  are  solved  simultane¬ 
ously,  and  then  the  constituent  concentration  is  solved.  This  method  improves 
solution  efficiency  dramatically  over  the  simultaneous  solution  of  all  four 
equations  and  unknowns. 

Hence,  the  solution  of  a  system  of  four  equations  and  four  unknowns 
becomes  the  solution  of  a  system  of  three  equations  (1,2,  and  8)  and  three 
unknowns  (m,  v,  and  h),  followed  by  the  solution  of  one  equation  (5)  and  one 
unknown  (5). 


Geometric  transform 

In  order  to  use  a  fixed  geometry  to  model  a  system  with  a  time  varying 
vertical  dimension  (depth)  it  is  convenient  to  use  a  geometric  transformation  to 
map  the  system  to  a  fixed  geometry. 
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Time  varying  system 


Fixed  grid  system 


The  transformation  is  based  on  the  following  relation: 

h  (11) 
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+  zb 


(12) 


Hence: 


U(x,y,z)  =  u 
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Y\ 
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l! 
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)) 

(13) 


After  completing  the  transformation  of  the  terms  and  simplifying,  we  arrive 
at  the  following  transformed  equations: 
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The  Momentum  equations 
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Volume  continuity 


=  0(15) 
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Advection-diffusion  equation 
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where 
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2-D  vertically  averaged  equations 

If  w,  v,  and  s  are  assumed  constant  with  respect  to  elevation  (z),  the  3-D 
equations  can  be  integrated  over  depth  to  yield  2-D  vertically  averaged 
equations.  For  example,  the  X-momentum  equation  reduces  to  the  following: 
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Similarly,  the  continuity  equation  reduces  to: 


du  Bv 
dx  dy 
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And  the  advection-diffusion  equation  reduces  to: 


h(b  -a)^-  +  h(b-  a)u^-  +  h(b  -  a)v^- 
dt  dx  dy 
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2-D  laterally  averaged  equations 

Lateral  averaging  eliminates  the  momentum  equation  in  the  direction  normal 
to  the  dominant  flow  direction.  The  equations  are  integrated  across  the  width  of 
the  channel.  This  operation  requires  that  the  channel  width  c  is  specified.  For 
the  purposes  of  TABS-MDS,  the  channel  width  in  laterally  averaged  elements  is 
constrained  such  that  it  is  constant  with  respect  to  depth,  but  can  vary  with 
respect  to  x  and  y  (i.e.,  along  the  channel  length).  For  example,  the  X- 
momentum  equation  reduces  to  the  following. 
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Similarly,  the  continuity  equation  reduces  to: 
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And  the  advection-diffusion  equation  reduces  to: 
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1-D  equations 

Under  this  approximation  both  vertical  and  lateral  integration  are  applied. 
Hence,  the  form  of  the  cross  section  must  be  defined.  In  TABS-MDS,  the  cross 
section  is  assumed  trapezoidal,  with  allowance  made  for  off-channel  storage.  For 
example,  the  X-momentum  equation  reduces  to  the  following: 
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A— +  Au¬ 
di  ox 
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Similarly,  the  continuity  equation  reduces  to: 


A 


rdu~]  ,  dA  d(A  +  A0C) 
—  +u—  + - r - 
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And  the  advection  diffusion  equation  reduces  to: 


\{A  +  Aoc)dl+Adx'AdX 
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Xdx 
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where 

A  =  main  channel  cross-sectional  area 
Aqc  =  off-channel  storage  cross-sectional  area 
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Finite  element  formulation 


To  generate  the  finite  element  equations,  integrate  each  of  the  equations  over 
the  element  volume  (for  3-D),  area  (for  2-D),  or  length  (for  1-D),  remembering  to 
include  the  weight  function  in  the  integration  (which,  for  the  Galerkin  method,  is 
the  same  as  the  basis  function). 

In  addition,  recast  the  higher-order  terms  using  integration  by  parts.  This 
causes  the  boundary  terms  to  drop  out  of  the  equations.  For  example,  take  the 
following  pressure  term,  multiplied  through  by  a  weight  function  N : 


N- 


pgh  dh 
(b  -  a)  dx 


(30) 


This  can  be  rewritten  as: 
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pg 
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(31) 


Then  ,  it  can  be  integrated  by  parts: 


N 


P  g  dlf_ 
2 (b-a)  dx 


d_ 

dx 


(n  ^  ) 

dN 

f  Pgh2  ) 

2  (b-a)) 

dx 

^2  (b-a)) 

-N 


gb1  dp 
2  (b-a)dx 


-N 


pgh2  da 
2(6 -a)2  ~dx 


(32) 


Note  that  the  first  term  on  the  right-hand  side  of  the  equation  can  be  evaluated  as 
an  area  integral  via  the  Gauss  Divergence  Theorem.  Hence,  it  becomes  a 
boundary  term. 


Time  derivative  solution  method 

The  time  derivative  is  approximated  with  a  simple,  fully  implicit  finite 
difference  formulation,  i.e., 


ap,  _(P,-(U) 

dt  At 
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where 


p,  =  any  of  the  unknown  variables  at  time  t 
At  =  the  time-step 


Newton-Rhapson  implementation 

Once  the  finite  element  equations  are  built,  they  are  solved  using  the 
Newton-Rhapson  iterative  method.  In  order  to  do  this,  partial  derivatives  with 
respect  to  each  of  the  unknown  variables  must  be  derived  for  each  system 
equation.  These  derivatives  compose  the  stiffness  matrix  and  are  used  to  dnve 
the  residual  (i.e.,  the  integral  of  each  equation  across  an  element)  to  0. 
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Expressions  for  Applied  Loads  and  Turbulent 
Mixing 


Bed  shear  stress 

The  bed  shear  stress  is  given  by  a  modified  form  of  Manning’s  Equation,  as 
given  by  Christensen  (1970).  Any  of  3  expressions  can  be  used,  depending  on 

the  instantaneous  value  of  the  depth/roughness  height  ratio  {  —} .  The 

expressions  are  as  follows  (given  for  the  X-direction  only): 
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All 


where 


tx  =  bed  shear  in  the  x-direction 
k  =  roughness  height 
d  =  local  depth 
v  =  local  velocity 
g  =  gravitational  constant 
p  =  density  of  water 

k  is  found  as  a  function  of  Manning’s  n  from  the  following  expression: 


k  = 


^  1.486  j 


(38) 


Wind  stress 

The  wind  stress  is  given  by  the  following  expression  (given  for  the  x- 
direction  only): 

T  =o  C  V2  cos0  (39) 

wr  ra  w  w  w  v  / 


where 

xwx  -  wind  stress  in  x-direction 
pa  =  density  of  air 
Vw  =  wind  velocity 

0„,  =  direction  from  which  the  wind  is  blowing,  measured  counterclockwise 
from  positive  x-axis. 

Cw  =  wind  stress  coefficient 
The  wind  stress  coefficient  is  given  by  Wu  (1980). 


0.8  + 0.065  xFm, 
1000 


(40) 
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Horizontal  turbulent  mixing  and  diffusion 

Horizontal  turbulent  mixing  can  be  specified  directly,  or  it  can  be  controlled 
by  the  method  of  Smagorinsky  (1963).  A  description  of  this  method  follows. 

The  Smagorinsky  method  of  describing  horizontal  eddy  viscosities  and 
diffusion  coefficients  is  a  “tensorially  invariant  generalization  of  the  mixing 
length  type  representation”  (Speziale  1998).  The  Smagorinsky  description  of  the 
turbulent  mixing  terms  in  the  Navier-Stokes  equations  are  given  as  follows.  For 
the  x-momentum  equation 
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For  the  y  momentum  equation 
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lc  =  Smagorinsky  coefficient,  usually  given  a  value  ranging  from  approxi¬ 
mately  0.005  for  rivers  to  0.05  for  estuaries  and  lakes  (Speziale  1998; 
Thomas  and  Williams  1995) 

A  =  surface  area  of  the  element 


The  Smagorinsky  description  of  the  turbulent  diffusion  terms  in  the  advection- 
diffusion  equation  is  given  as  follows: 


dx  V 
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dx 
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In  order  to  promote  numerical  stability,  TABS-MDS  provides  a  means  of 
establishing  minimum  values  of  turbulent  mixing  and  turbulent  diffusion.  These 
values  are  used  in  place  of  the  Smagorinsky  term  (5)  when  they  are  found  to 
exceed  the  value  of  that  term.  The  minimum  turbulent  mixing  value  is  given  by 
the  following  equation: 
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Se  min  —  TBMINF  x  pa-jA 


(45) 


The  minimum  turbulent  diffusion  value  is  given  by  the  following  equation: 


Sd  min  =  TBMINFSxocVl’ 


(46) 


where 


TBMINF  =  minimum  turbulent  mixing  factor  (default  =  1 .0) 
TBMINFS  =  minimum  diffusion  factor  (default  =  1 .0) 


a  =  a  coefficient,  given  as  5.00xl0'3  ft/sec  or  1.52xl0'3  m/s, 
depending  on  the  unit  system  being  used  in  the  simulation. 

This  value  is  an  arbitrary  estimate  of  the  minimum  turbulent 
mixing  needed  to  ensure  model  stability.  It  equals  the  value  of 
eddy  viscosity/diffusion  which  corresponds  to  a  Peclet  number 
of  40  and  a  velocity  magnitude  of  0.2  ft/sec. 


Also,  if  |  V  |  <  TBMINF  x  Vmm,  S/.:u-lU1  is  applied,  regardless  of  the  turbulent 
mixing  as  given  by  the  Smagorinsky  calculation.  This  is  done  to  inhibit 
numerical  instability  in  areas  with  both  extremely  small  velocities  and  high 
velocity  gradients. 


Vertical  turbulent  mixing  and  diffusion 

Vertical  turbulent  mixing  and  diffusion  are  given  by  the  method  of  Mellor- 
Yamada  (1982)  with  a  modification  according  to  Hendersen-Sellers  (1984). 

The  Mellor-Yamada  expressions  for  vertical  eddy  viscosity  and  diffusion  are 
given  as  follows: 


Ea=Eyz=pSmlmq 
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(50) 
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Sm  =  0.393 
Sh  =  0.494 
bi  =  16.6 


The  Henderson-Sellers  adjustment  is  a  factor  that  accounts  for  the  dampening 
affect  on  turbulence  induced  by  stable  stratification.  The  factor  is  expressed  in 
terms  of  the  Richardson  number: 


/?,  =- 


-g(3p/3z) 
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du 
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For  vertical  diffusion  of  momentum  (i.e.,  vertical  eddy  viscosity)  the  expression 
is  given  as  follows: 


(1  +  0.74/?,.) 


(52) 


where  Ez  is  the  vertical  eddy  viscosity,  and  E:o  is  the  vertical  eddy  viscosity 
assuming  no  stratification  influence  on  the  turbulence  (i.e.,  the  value  taken  from 
Mellor-Yamada). 

For  vertical  diffusion  of  salinity  (i.e.,  vertical  diffusion  coefficient)  the 
expression  is  given  as  follows: 


D.  = 


D 


(1  +  37/?/) 


(53) 


where  Dz  is  the  vertical  diffusion  coefficient,  and  Dzo  is  the  vertical  diffusion 
coefficient  assuming  no  stratification  influence  on  the  turbulence  (i.e.,  the  value 
taken  from  Mellor-Yamada). 
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